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Abstract 


Most  objects  tracked  in  space  follow  a  regular  Keplerian  orbit;  unfortunately, 
non-Keplerian  objects  such  as  maneuvering  satellites,  tethered  systems,  and  thrusting 
ballistic  missiles  are  becoming  more  common.  It  is  important  to  be  able  to  distinguish 
between  Keplerian  and  non-Keplerian  objects  due  to  the  potential  risk  of  a  tethered 
satellite  being  mistaken  for  an  object  on  re-entry.  This  research  focused  on  creating  a 
computer  model  that  can  detect  the  non-gravitational  acceleration  present  in  non- 
Keplerian  orbits.  A  3rd  order  Taylor  series  expansion  was  used  to  model  the  dynamics 
and  to  produce  simulated  radar  data.  Linear  least  squares  estimation  was  used  to  estimate 
the  initial  state  of  a  space  object  with  a  state  vector  composed  of  position,  velocity, 
acceleration,  and  its  first  derivative.  Monte  Carlo  analysis  was  used  to  verify  that  the 
estimator  was  unbiased  and  representative  of  the  uncertainty  in  the  data.  The  Monte 
Carlo  method  detected  non-gravitational  acceleration  as  small  as  1.12  cm/s  ;  however,  a 
subsequent  approach  that  analyzed  the  data  sets  individually  only  detected  acceleration  as 
small  as  10.63  cm/s".  At  smaller  magnitudes,  the  estimator  was  able  to  detect  the 
presence  of  non-gravitational  acceleration,  but  was  ultimately  unable  to  estimate  the  true 
value  with  statistical  accuracy. 
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RADAR  ORBIT  ANALYSIS  TOOL 
USING  LEAST  SQUARES  ESTIMATOR 


I.  Introduction 


Background 

Radars  have  been  used  to  track  satellites  ever  since  the  Soviet  Union  successfully 
launched  Sputnik  in  October  1957.  These  radars,  along  with  optical  sensors,  eventually 
grew  into  what  is  commonly  known  today  as  the  Space  Surveillance  Network  (SSN). 

The  SSN  has  25  sensor  sites  located  world-wide  with  the  capability  of  tracking  space 
objects  from  low  earth  orbit  (LEO)  all  the  way  to  geosynchronous  earth  orbit  (GEO). 
There  are  currently  over  10,000  objects  being  tracked  today;  however,  the  number  of 
objects  increases  each  year  due  to  more  launches  and  the  increased  capability/desire  to 
track  smaller  and  smaller  objects.  These  objects  include  satellites,  rocket  debris,  bolts, 
and  sometimes  even  an  astronaut’s  glove. 

The  various  sensors  used  within  the  SSN  offer  their  own  capabilities  and 
challenges.  Optical  sensors  use  a  Charge  Coupled  Device  (CCD)  to  take  pictures  of  the 
satellite  streaking  across  the  sky.  These  images  are  cross-checked  with  star  charts  to 
determine  the  satellite’s  topocentric  right  ascension  and  declination  (Vallado,  2001:243). 
Optical  sensors,  unfortunately,  are  not  practical  during  the  day  nor  do  they  work  very 
well  when  it  is  cloudy  or  overcast. 

Radars,  however,  can  be  used  day  or  night  and  in  cloudy  conditions.  Radars  use 
both  a  transmitter  to  send  electromagnetic  radiation  and  a  receiver  to  obtain  the  EM 
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energy  reflected  back  from  the  detected  object.  This  method  obtains  the  topocentric 
azimuth  and  elevation  of  the  space  object.  The  SSN  divides  radars  into  three  categories: 
tracking,  detecting  and  phased  array  (interferometers).  Interferometers  use  a  collection  of 
transmitters  and  receivers  to  detect  objects.  The  transmitters  create  a  fan  of  energy  that  is 
reflected  back  to  the  receivers  when  an  object  passes  through.  Varying  the  phase  creates 
the  ability  to  obtain  direction  cosines  that  can  be  manipulated  to  obtain  elevation  and 
azimuth  data  (Vallado,  2001:241). 

Radars  used  for  tracking  are  typically  more  accurate  than  those  used  for  initial 
detection  (Thomas,  1967:75-86);  however,  all  systems  contain  a  certain  level  of  bias  and 
noise.  Biases  in  range,  azimuth,  and  elevation  are  taken  into  account  to  correct  the  raw 
sensor  data.  Noise  statistics  are  very  useful  in  helping  to  accurately  weigh  the  validity  of 
the  data.  This  is  important  when  the  data  is  to  be  run  through  estimators  or  filters. 

In  the  50  years  since  Sputnik,  the  ability  to  track  space  objects  has  become 
increasingly  more  complex  due  to  technological  advancements  both  foreign  and 
domestic.  Satellites  can  follow  more  than  just  a  regular  Keplerian  orbit.  For  a  given 
radar  track,  the  space  object  may  be  maneuvering,  part  of  a  tethered  system,  or  even  a 
thrusting  ballistic  missile. 

Tethered  systems  are  of  key  interest  these  days.  The  force  created  by  a  tether 
changes  the  dynamics  of  the  end  masses  and  may  result  in  the  satellite  being  mistaken  for 
an  object  on  a  re-entry  trajectory  (Asher  and  others,  1988:5 14).  If  the  end  mass  is  in  a 
lower  orbit  than  the  center-of-mass,  its  velocity  will  be  smaller  than  what  should  be 
expected  from  a  single  satellite  in  a  Keplerian  orbit;  similarly,  an  end  mass  in  a  higher 
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orbit  than  the  center-of-mass  will  have  a  higher  velocity  than  it  should  (Cicci  and  others, 
2002:340). 

Problem  Statement 

In  order  to  detennine  if  an  object  is  following  a  regular  Keplerian  orbit  or  has 
some  other  motion,  it  is  necessary  to  determine  how  the  dynamics  between  the  two 
groups  differ.  One  solution  is  to  develop  a  model  that  can  detect  non-gravitational 
acceleration.  This  model  must  be  easily  implemented  and  have  the  ability  to  produce 
results  relatively  quickly.  Such  a  model  can  be  produced  using  Galileo’s  projectile 
motion  dynamics  and  a  linear  least  squares  estimator.  A  linear  least  squares  estimator  is 
ideal  because  there  is  no  need  for  an  initial  guess  or  iteration,  so  results  are  produced  in  a 
timely  manner. 

Research  Objectives 

The  primary  objective  of  this  research  was  to  create  a  statistically  accurate 
computer  model  capable  of  determining  if  a  space  object  has  non-gravitational 
acceleration.  A  truth  model  and  least  squares  estimator  must  be  produced  in  order  to 
create  and  test  the  computer  model.  This  model  is  merely  a  filter  to  be  used  at  radar  sites 
to  get  an  initial  idea  if  a  space  object  is  following  a  non-Keplerian  orbit. 

Research  Focus 

This  research  focused  on  using  the  projectile  equations  of  motion  with  a  linear 
least  squares  estimator  to  solve  the  estimation  problem  instead  of  a  sequential  method 
such  as  either  a  Bayes  filter  or  a  Kalman  filter.  The  state  vector  will  include  not  only 
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position  (r)  and  velocity  (r ),  but  also  acceleration  (r )  and  its  first  derivative  (r ).  The 
inclusion  of  acceleration  in  the  state  vector  was  used  to  verify  if  non-gravitational 
acceleration  was  present.  The  covariance  must  also  statistically  validate  that  the 
acceleration  exists. 

Investigative  Questions 

Several  questions  were  addressed  during  this  research.  The  most  important 
question  was,  “What  magnitude  of  acceleration  can  the  model  detect?”  Tethered  systems 
and  maneuvering  space  objects  can  have  non-gravitational  accelerations  that  are  quite 
small  (i.e.  cm/s",  pm/s“).  The  computer  model  will  be  more  versatile  if  it  can  detect  these 
smaller  accelerations.  Another  question  of  interest  was,  “Are  the  dynamics  accurate 
enough?”  Galileo’s  equations  for  projectile  motion  are  very  general.  A  Taylor  series 
expansion  of  Galileo’s  equations  of  motion  may  be  utilized  for  better  accuracy.  The 
question  of  whether  or  not  the  geopotential  is  accurate  enough  with  just  h  and  two-body 
tenns  was  also  addressed. 

Methodology 

Solving  the  estimation  problem  required  dividing  the  process  into  four  stages. 
First,  the  raw  radar  data  was  converted  to  a  more  usable  format.  Second,  a  truth  model 
that  can  simulate  various  aspects  of  real-world  data  was  developed.  Third,  an  estimator 
that  outputs  r ,  r,  r,  and  r  was  created.  The  estimator  must  also  produce  a  covariance 
matrix  to  verily  that  the  state  estimate  was  accurate.  Last,  test  case  scenarios  were  run 
through  the  truth  model  and  estimator  to  ensure  proper  function  and  accuracy. 
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Assumptions  and  Limitations 

A  number  of  assumptions  and  limitations  were  taken  into  account  and  addressed 
during  this  research.  The  equations  of  motion  used  for  the  truth  model  were  produced 
using  Galileo’s  projectile  trajectory  with  a  Taylor  series  expansion.  Two-body  equations 
of  motion  are  more  accurate  in  modeling  space  object  motion;  however,  those  equations 
are  nonlinear  and  would  make  the  estimation  problem  more  complex.  For  these 
equations  of  motion,  the  value  for  gravity  includes  more  than  the  standard  9.81  m/s“. 
Gravity  was  instead  modeled  as  the  gradient  of  the  geopotential  taking  into  account  both 
two-body  and  J2  effects.  Air  drag  was  not  modeled.  The  omission  of  air  drag  works  well 
for  a  short  radar  track,  but  would  limit  the  accuracy  of  the  model  if  used  for  a  much 
longer  simulation.  For  this  research,  the  local  gravity  vector  was  calculated  only  once, 
using  the  initial  position  vector.  For  such  a  short  track  of  data,  it  was  assumed  that  the 
gravity  vector  would  not  change  dramatically. 

This  model  is  limited  by  radar  capabilities  and  how  fast  the  space  object  is 
accelerating.  Assuming  the  radar  has  a  range  error  of  100  meters,  given  a  five-minute 
radar  track,  the  object  must  be  accelerating  faster  than  0.667  cm/s“  in  order  for  the  model 
to  statistically  say  there  is  an  extra  acceleration  component  present.  This  model  is, 
therefore,  unable  to  detect  extremely  small  accelerations.  It  was  also  assumed  that  the 
sensor’s  total  instrument  error  was  composed  of  only  little  error  sources,  thereby 
allowing  the  concept  of  Gaussian  Distribution  to  be  taken  into  account  during  estimation. 
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Implications 

This  computer  model  is  a  simple  tool  that  radar  sites  run  their  data  through  to 
quickly  get  an  idea  if  a  space  object  is  following  a  non-Keplerian  orbit.  It  is  merely  a 
filter.  It  can  determine  if  non-gravitational  acceleration  is  present;  however,  it  cannot 
distinguish  between  a  tethered  system  and  a  maneuvering  satellite.  A  more  complex 
model  would  need  to  be  utilized  to  provide  specific  object  identification  and  orbit 
propagation. 

Preview 

Chapter  II  touches  on  some  of  the  key  research  areas  studied  while  solving  the 
problem  at  hand,  the  most  important  of  which  being  least  squares  and  how  it  all  works. 
The  importance  of  Galileo’s  projectile  trajectory  is  also  examined.  A  summary  of 
previous  attempts  to  model  space  objects  with  additional  acceleration  components  is 
included.  Chapter  III  discusses  the  specific  processes  used  to  create  the  conversion  from 
raw  radar  data,  the  truth  model,  and  the  estimator;  as  well  as  how  the  routines  are 
executed  and  validated.  Chapter  IV  goes  over  the  different  test  case  scenarios,  as  well  as 
an  analysis  of  the  data.  Chapter  V  concludes  with  the  end  result  and  recommendations 
for  future  research. 
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II.  Literature  Review 


Chapter  Overview 

The  purpose  of  this  chapter  is  to  explain  the  theory  behind  the  key  concepts  used 
in  this  research.  Section  2.2  covers  some  of  the  previous  research  efforts  related  to  space 
objects  with  non-Keplerian  orbits.  Section  2.3  explores  Galileo’s  insight  in  projectile 
motion  and  how  it  can  be  adapted  to  modeling  satellite  motion.  Section  2.4  explains  the 
fundamental  assumptions  required  to  make  least  squares  work,  as  well  as  how  to  obtain 
the  estimated  state  and  covariance  for  a  set  of  data. 

Non-Keplerian  Orbits 

Orbit  determination  of  space  objects  has  been  a  topic  of  interest  ever  since 
astronomers  first  tried  tracking  the  planets  and  moons  in  the  solar  system.  Many 
scientists  and  astronomers,  including  Galileo,  Brahe,  Kepler,  and  Newton  provided 
insight  into  how  to  track  these  objects.  It  was  Newton’s  work  combined  with  Kepler’s 
Planetary  Laws  that  created  the  two-body  equations  of  motion: 

n  fir 

r= - r  (1 

r 

Equation  (1)  is  fundamental  in  understanding  the  dynamics  of  orbiting  objects.  Many 
different  orbit  detennination  methods  use  the  two-body  equations  of  motion  as  the 
foundation. 

Also  noteworthy  was  the  development  of  Kepler’s  Problem.  Given  an  initial 
position  and  velocity  of  a  space  object,  plus  a  time  span  (t  —  t0),  the  position  and  velocity 
at  the  future  time  could  be  detennined.  This  method  revolutionized  the  orbit 
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determination  process.  The  propagation  of  objects  following  a  Keplerian  orbit  was  now 
obtainable.  A  thorough  discussion  on  Kepler’s  Problem  can  be  found  in  Fundamentals  of 
Astrodynamics  (Bate  and  others,  1971:177-203)  ox  Fundamentals  of  Astro  dynamics  and 
Applications  (Vallado,  2001:87-103). 

Kepler’s  work  has  been  exceedingly  useful  in  tracking  satellites  orbiting  Earth. 

The  development  of  tethered  systems,  however,  has  presented  a  problem:  the  end  masses 
of  tethered  systems  do  not  follow  a  regular  Keplerian  orbit.  The  tether  creates  an 
additional  acceleration  in  the  radial  and  tangential  directions.  The  observed  ‘daughter’ 
satellite  (m)  and  the  unobserved  ‘parent’  satellite  (mp)  with  their  force  components  are 
presented  below  in  Figure  1. 


Figure  1.  Tether  Force  Components  (Cicci  and  others,  2001a:314) 


As  a  result  of  this  new  development,  a  plethora  of  orbit  determination  methods  for 
tethered  systems  have  been  researched  over  the  past  decade.  Most  techniques  used  to 
model  tethered  orbits  assume  the  following  equations  of  motion 


"  -ur  F 
r  -  — 3-  +  — 
r  m 


(2) 
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(3) 


,  -pr  ^  F 

rp=^^  +  ~ 
rP  mP 

Many  attempts  to  model  tethered  systems  have  focused  on  estimating  the 
additional  acceleration  components  seen  in  Figure  1 .  Thus,  the  state  vector  is  commonly 
written  as 


x 


y 

z 


X  = 


x 

y 


z 


(4) 


Engineers  have  used  least  squares  (Cicci  and  others,  2001a:309-326)  and  ridge-type 
estimators  (Cicci  and  others,  2001b:297-316)  to  estimate  the  initial  state  of  tethered 
systems.  More  in-depth  techniques  have  been  used  to  find  libration  angle,  libration  rate, 
and  trajectory  prediction  (Cicci  and  others,  2002:340).  These  techniques  use  more 
complex  dynamics  and  require  longer  arcs  of  observation  data.  Since  the  equations  of 
motion  listed  above  are  non-linear,  the  estimation  process  requires  an  initial  guess  for  the 
state. 

It  is  not  just  tethered  systems  that  contain  additional  acceleration  components. 
Maneuvering  space  objects  produce  additional  acceleration  that  creates  a  non-Keplerian 
orbit.  Satellites  use  thrust  for  various  reasons,  such  as:  station  keeping,  rendezvous,  and 
orbital  transfers.  Thrust  is  a  non- instantaneous  force  that  can  be  written  as  F=ma,  using 
Newton’s  Law;  therefore,  when  thrust  is  present,  it  is  an  additional  acceleration  that 
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ultimately  causes  a  non-Keplerian  orbit.  Since  this  research  primarily  focused  on  the 
orbit  determination  of  any  non-Keplerian  object,  not  just  tethered  systems,  the  equations 
of  motion  and  estimation  techniques  mentioned  above  will  be  replaced  with  a  more 
general  model. 


Galileo’s  Projectile  Trajectory 

In  1638,  Galileo  Galilei  helped  to  redefine  man’s  understanding  of  motion. 
Galileo’s  book,  Dialogues  Concerning  Two  New  Sciences  (Galilei,  1914),  covered  four 
days  of  experimentation  and  contemplation  that  ultimately  created  the  foundation  for 
projectile  equations  of  motion.  Galileo  conveyed  the  following  insights  (Hahn, 
2002:341): 

1 .  All  bodies  falling  in  a  vacuum  do  so  with  the  same  constant 
acceleration.  For  a  body  falling  from  rest,  the  speed  is  proportional  to 
the  elapsed  time.  This  is  so  both  in  the  situation  of  free  fall  and  for 
balls  rolling  on  an  inclined  plane. 

2.  The  law  of  fall,  namely,  that  the  distance  covered  by  a  body  moving 
from  rest  is  proportional  to  the  square  of  the  time  of  the  motion. 

3.  The  trajectory  of  a  projectile  has  a  parabolic  shape. 

It  is  not  until  later,  with  the  development  of  calculus,  that  Galileo’s  insights  were  put  into 
the  familiar  fonn 


.y 

II 

1 

Qrq 

(5) 

y  =  vyt-\gt2 

(6) 

Vv  =  V*0 

(7) 

X  =  Vj 

(8) 
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It  was  Galileo’s  work  concerning  the  motion  of  projectiles  that  produced  the  basic 
equations  of  motion  used  for  this  research.  Galileo  studied  the  properties  of  both  uniform 
motion  and  naturally  accelerated  motion  and  proposed  that  a  projectile  is  a  combination 
of  the  two.  The  horizontal  motion  is  produced  by  the  launching  mechanism,  whereas  the 
vertical  motion  is  due  to  gravity.  A  fired  cannonball  is  a  prime  example  of  projectile 
motion. 

Satellites  and  other  space  objects  can  also  be  modeled  using  projectile  motion 
(although  a  more  accurate  model  would  use  two-body  equations  of  motion).  Figure  2 
helps  illustrate  the  basics  of  orbital  motion. 


Figure  2.  Trajectories  (Sellers,  undated:31) 

Imagine  a  person  standing  on  one  side  of  the  Earth  throwing  a  baseball.  The 
faster  the  ball  is  thrown,  the  further  it  travels.  Earth,  however,  is  not  flat:  it  is  round. 
Therefore,  as  the  ball  is  flying  through  the  air,  the  Earth  is  actually  curving  away  from  the 
ball  at  a  rate  of  5  meters  for  every  8  kilometers  traveled.  At  the  right  velocity,  the  object 
is  falling  slower  than  the  rate  at  which  the  Earth  is  curving  away.  At  this  point,  the  object 
has  reached  ‘freefall,’  also  known  as  orbit.  An  object  must  be  moving  at  7.9  km/s 
(ignoring  air  drag)  to  be  considered  in  orbit  (Sellers,  2000: 106).  With  this  example  in 
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mind,  it  is  easy  to  see  how  Galileo’s  first  insight  listed  above  is  applicable  to  space 
objects. 


In  order  to  find  the  “constant  acceleration”  exerted  on  the  space  object,  Newton’s 

d  ,  -  F 

2  Law  is  required.  Rearranging  Newton’s  2  Law  yields  a  = —  .  In  space,  where 

m 


weight  (W)  =  force  (F),  Newton’s  law  looks  like 


r_F_  W_ 
m  m 


Substituting  mg  for  W  and  simplifying  yields 


a  =  g 


(9) 


(10) 


The  above  discussion  has  just  shown  that  space  objects  are  subject  to  constant 
acceleration  in  the  form  of  gravity  as  long  as  gravity  is  assumed  to  be  constant.  It  seems 
only  reasonable,  therefore,  to  be  able  to  model  objects  in  LEO  using  projectile  motion, 
with  the  understanding  that  the  model  is  very  crude  and  only  ‘good  enough’  for  a  short 
period  of  time.  Below  are  Galileo’s  equations  in  vector  form. 


r  =  r0+v0At  +  ^gAt2 

(ID 

v  =  v0  +  gAt 

(12) 

a  =  g 

(13) 

A  quick  glance  at  the  above  equations  verifies  that  they  are  actually  a  2nd  order  Taylor 
series  approximation.  If  the  2nd  order  does  not  yield  accurate  results,  then  an  expansion 
to  the  3rd  order  or  even  4th  order  may  prove  better.  The  equations  of  motion  for  a  3rd 
order  and  4th  order  Taylor  series  approximation  are  listed  in  Appendix  A. 
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Least  Squares 

In  1801,  German  mathematician  Carl  Friedrich  Gauss  discovered  a  technique  that 
would  revolutionize  the  methods  of  orbit  determination.  Gauss’  discovery,  known  as 
Least  Squares,  relies  on  several  principle  assumptions:  the  dynamics  contain  no  error,  the 
instruments  do  contain  error,  Gaussian  distribution,  and  Principle  of  Maximum 
Likelihood  (Hall,  1994:7).  In  order  for  the  dynamics  to  contain  no  error,  it  is  important 
to  model  the  object  of  interest  with  a  relevant  dynamics  system.  The  Central  Limit 
Theorem  addresses  the  instrument  error  and  Gaussian  distribution.  The  Principle  of 
Maximum  Likelihood  is  the  final  assumption  required  in  order  to  produce  the  estimate  of 
a  state  and  its  covariance.  These  assumptions  are  addressed  in  the  following  sub-sections. 

Central  Limit  Theorem. 

The  Central  Limit  Theorem  states  that  if  an  instrument  has  many  little  errors,  then 
no  matter  how  the  little  errors  are  distributed,  the  overall  error  can  be  described  using  a 
Gaussian  function  (Figure  3). 


Figure  3.  Gaussian  Distribution  Function  (Zaninetti,  2002) 

The  Gaussian  function  is  centered  about  the  true  value.  The  width  of  the  curve  is 
described  by  the  standard  deviation  ( o ).  A  smaller  o  means  a  narrower  curve.  The 
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probability  that  the  answer  is  within  ±  1  a  is  68%.  Within  ±2  <7  yields  a  probability  of 
~95  %  while  ±  3  o  is  99.7%.  The  Gaussian  Distribution  can  be  written  as 


2  c2 


(14) 


The  Gaussian  Distribution  is  an  important  concept  utilized  by  statisticians,  scientists,  and 
engineers. 

Expectation  Operator. 

The  expectation  operator  is  a  linear  operator  that  facilitates  the  estimation  process. 
The  expectation  operator  can  be  written  as 


E{-)=](-)f(X)dX 


(15) 


The  term  f{X)  is  the  probability  density  function  and  (-)  is  the  variable  of  interest. 

After  taking  Equation  (14)  and  inputting  it  into  Equation  (15),  the  expectation  operator 
looks  like 

E(~)  =  \  (->— ^<2fr2°  dX  06) 

crv2/r 

Equation  (16)  yields  several  cases  of  interest: 

E{X)  =  X o  (17) 

E(X-X0)  =  0  (18) 

E((X-X0)2)  =  E(e2)  =  c J2  (19) 

These  results  imply  that:  “the  expected  value  of  a  measurement  is  the  true  value” 
(Equation  (17)),  “the  average  error  in  the  measurement  is  zero”  (Equation  (18)),  and  “the 
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average  squared  error  is  a1  ’’(Equation  (19))  (Wiesel,  2003b:  19).  The  result  of  Equation 
(19)  plays  a  key  role  in  least  squares  estimation. 

Principle  of  Maximum  Likelihood. 

Suppose  that  N  measurements  are  taken  of  an  object.  As  long  as  each 
measurement  is  independent,  the  probability  of  obtaining  the  data  set  is  the  product  of  the 
individual  probability: 


P(X],X2...XN)  = 


(2k) 


N/2 


JL  1 

riT 


M 


e  i=1 


,(Xt-X0? 
I  2  of 


(20) 


Unfortunately,  no  matter  how  well  the  object  of  interest  is  measured,  xo  cannot  be 
obtained.  Therefore,  a  different  approach  must  be  attempted. 

The  Principle  of  Maximum  Likelihood  is  where  the  estimate  X  is  defined  as, 
“the  value  of  X0  which  maximizes  the  probability  of  having  obtained  the  actual  data  set’' 
(Wiesel,  2003b:20).  Equation  (20)  now  looks  like 


P(XVX2...XN)  = 


(2k) 


N/2 


JL  1 

fit 


WV-AT 
f  2  c} 


(21) 


Subsequently,  the  true  error  (e  =  X  —  X0)  has  been  exchanged  for 


Equation  (21)  is  maximized  when  the  tenn  within  the  exponential 


d_f(Xi-X)2_Q 
dX  2<t2 


a  residual  (  r  =  X  -  X  ). 
is  minimized: 


(22) 


which  yields 


i=i  °i 


(23) 


and  simplifies  to 
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(24) 


i= 1 

where 

1 


(25) 

The  above  process  of  minimizing  the  exponential  is  how  the  Method  of  Least 
Squares  acquired  its  name.  As  stated  by  Weld  (1916:59),  “the  most  probable  value  of  a 
measured  quantity  that  can  be  deduced  from  a  series  of  direct  observations,  made  with 
equal  care  and  skill,  is  that  for  which  the  sum  of  the  squares  of  the  residuals  is  a 
minimum.” 

As  with  any  answer  that  is  detennined,  it  is  important  to  know  how  accurate  the 
answer  really  is;  therefore,  the  variance  of  X  must  be  found: 

<72=£((X-X0)2)  (26) 

After  inserting  Equation  (24)  into  the  above  equation,  rearranging,  and  simplifying,  the 
result  is 


W,  =  ■ 


I 


°7 


al  =- 


N  1 

zi 


(27) 


Using  the  result  from  Equation  (27),  the  estimate  from  Equation  (24)  can  then  be 
simplified  to 


x=<?;Z 


X, 


(28) 
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Equations  (27)  and  (28)  are  the  two  most  important  equations  of  the  estimation  process. 
Given  data  Xi  and  its  standard  deviation  cr  ,  the  estimate  of  the  true  value  and  its 
standard  deviation  is  now  obtainable. 

Multi-Dimensional  Probability. 

Thus  far,  the  Gaussian  function  has  been  written  as  the  one  variable  case.  In  orbit 
determination,  the  Gaussian  function  will  be  multi-dimensional: 


where 


f{X) 


N/ 

{7.71)  72 


1/  T(*-*o)  P 

1  72  eK 


(29) 


(30) 


The  covariance  matrix  (P)  is  a  positive  semi-definite  matrix.  The  diagonal  terms  are  the 
(J2  quantities  and  are  called  the  variances.  It  is  the  square  root  of  the  variance  that  relates 
the  accuracy  of  the  estimate  of  the  state.  The  covariance  matrix  is  normally  defined  as 

P  =  E(eeT )  (31) 

The  above  equation  looks  remarkably  similar  to  Equation  (19)  from  the  one  variable 
problem.  This  is  why  the  diagonal  terms  of  the  covariance  matrix  are  defined  as  o2 . 

Linearized  Dynamics. 

Unlike  the  above  process  where  just  one  component  of  an  object  is  estimated, 
orbit  determination  is  usually  composed  of  several  different  components.  These 
components;  such  as  position,  velocity,  and  acceleration,  are  placed  into  a  state  vector 
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(x) .  Engineers  and  scientists  are  interested  in  how  the  state  of  a  space  object  changes 


with  time.  This  can  be  written  as 


dX 

dt 


g{X,t) 


(32) 


or 

X(t)  =  h(X(t0),t)  (33) 

Equation  (33)  shows  the  actual  solution  written  in  terms  of  the  initial  state  and  time.  The 
state  transition  matrix  ( O )  propagates  the  actual  state  as  a  function  of  time.  If  the 
dynamics  can  be  written  as  Equation  (33),  then  O  can  be  written  as 

O(t,t0)  =  VxWh(X(t0),t)  (34) 


or 


®(t,t  0) 


a x(t) 

dX(t0) 


(35) 


This  approach  greatly  simplifies  the  estimation  process. 

Linear  Least  Squares. 

Linear  least  squares  is  an  estimation  process  used  on  systems  where  the  dynamics 
can  be  written  as  linear  differential  equations  ( i.e .  the  function  and  its  derivatives  appear 
only  in  their  first  power  and  are  also  not  multiplied  with  each  other).  The  observations 
must  also  be  linear  in  order  to  utilize  linear  least  squares.  Orbit  determination  problems 
are  generally  nonlinear  because  they  use  two-body  equations  of  motion  to  model  the 
dynamics;  and  the  observations  are  in  the  fonn  of  range,  azimuth,  and  elevation,  which 
are  also  nonlinear. 
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The  state  of  a  linear  system  can  be  found  at  any  time  using  the  state  transition 


matrix: 


X(t)  =  O(t,t0)X(t0)  (36) 

As  long  as  the  observations  are  linear,  they  can  be  written  using  the  observation  relation : 

Zi(ti)  =  HiX{ti)  +  ei  (37) 

Substituting  Equation  (36)  into  the  above  equation,  and  solving  for  the  error  yields 

et=Zt  (O-H&h,  t0)X(t0)  (38) 

The  term  HiO(ti,t0 )  can  be  replaced  with  71  to  yield 

ei=Zi(ti)-TiX(t0)  (39) 

Given  that  true  error  can  never  be  determined,  the  above  equation  is  replaced  with  the 
residual: 


r=Z-TX(t0 ) 


(40) 


Before  proceeding,  the  following  shorthand  notation  is  used  to  simplify  the  least  squares 
equations: 


~TX~ 

"  77,0(6,0)  " 

1 2 

= 

772O(t2,0) 

Tn\ 

_HNQ(tN,0) 

(41) 


(42) 
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(43) 


Q  = 


Qi  o 
0  q2 


0 

0 


0  0  ...  qn 


where  Q  is  the  instrumental  covariance  matrix  and  its  inverse  is 


Q~'  = 


q r1  o 

o  Qx 

o  o 


o 

0 

q 


The  multi-dimensional  probability  function  resembles 

1  ,  ,-1/  \^(Z-TX(t0)fQ'(Z-TX(t0)) 

XX)  = - w\Q\/2e{2 

(2  7T)/2 


(44) 


(45) 


After  minimizing  the  exponential,  the  estimated  state  vector  at  t0  can  be  written  as 


X(0)  =  (TtQ~xT) TtQ aZ  (46) 

The  covariance  of  the  estimate,/5^  =  E(eye'Y ) ,  goes  through  several  lines  of  substitution 
and  simplification  to  arrive  in  the  fonn 

px  =  (TtQ-xT)-x  (47) 

The  state  of  a  system  and  its  covariance  can  now  be  estimated  given  a  batch  of  data  and 
the  instrumental  covariance. 


Summary 

The  ability  to  identify  and  track  space-based  objects  has  played  a  vital  role  in  the 
United  States’  space  superiority.  Scientists  and  engineers  have  utilized  the  concept  of 
Keplerian  motion  to  aid  in  tracking  and  propagating  space  objects;  however,  as 
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technology  improves  and  space  objects  become  more  complex  ( i.e .  increased 
maneuverability,  tethered  systems,  etc.),  the  space  objects  will  not  always  follow 
Keplerian  motion.  Many  new  methods  have  been  designed  to  track  non-Keplerian  orbits. 
Most  techniques  use  non-linear  equations  of  motion  for  the  dynamics  and  focus  on 
finding  the  additional  acceleration  components  created  by  either  the  tether  of  a  tethered 
satellite  system  or  the  thrust  from  a  maneuvering  space  object. 

Galileo’s  work  on  projectile  motion  lends  itself  to  linear  equations  of  motion  that 
can  be  adapted  to  modeling  objects  in  LEO.  These  equations  are  very  general  and  do  not 
take  into  account  atmospheric  drag;  also,  gravity  is  assumed  to  be  constant.  When 
modeling  orbits  with  these  equations  of  motion,  the  results  are  only  accurate  for  a  short 
arc  of  observation  data.  Linear  least  squares  can  be  combined  with  these  equations  of 
motion  to  estimate  the  state  of  the  system.  The  linear  least  squares  method  has  been  used 
for  hundreds  of  years  and  is  a  vital  part  of  estimation  theory. 
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ITT.  Methodology 


Chapter  Overview 

The  purpose  of  this  chapter  is  to  explain  the  processes  used  to  create  the 
MATLAB  computer  model.  Section  3.2  explores  how  to  transform  the  raw  radar  data 
into  position  vectors  for  use  in  the  estimator.  Section  3.3  explains  how  the  truth  model 
was  created  using  Galileo’s  projectile  motion  equations.  Section  3.4  explains  the  creation 
of  the  linear  least  squares  estimator,  which  produces  the  estimate  of  the  state  and  its 
covariance. 

Raw  Data 

Radar  data  is  presented  in  various  forms  depending  on  the  type  of  sensor  that 
obtains  it.  For  the  purposes  of  this  research,  it  was  assumed  that  the  radar  data  will 
consist  of  range  (p),  azimuth  (a),  and  elevation  (P)  values.  Range  is  measured  in 
kilometers  while  azimuth  and  elevation  are  measured  in  radians.  These  radar  values 
represent  the  space  object’s  position  in  the  radar  site’s  topocentric  frame,  where  the  axes 
are  labeled  South,  East,  and  Zenith  (SEZ).  Figure  4  represents  the  geometry  of  the 
problem.  Azimuth  is  measured  clockwise  from  the  north.  Elevation  is  measured 
upwards  from  the  site’s  horizon,  and  range  is  measured  from  the  radar  site  to  the  satellite 
or  space  object  being  observed. 
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Figure  4.  Observation  Geometry 


The  SEZ  frame  is  not  inertial;  it  rotates  with  the  Earth.  The  raw  data,  therefore, 
must  be  converted  from  the  SEZ  coordinate  frame  to  position  vectors  in  the  Earth 
Centered  Inertial  (IJK)  frame,  which  requires  several  steps.  Figure  5  helps  to  visualize 
the  problem. 


Figure  5.  Satellite  Position 


Rsite  represents  the  position  from  the  center  of  the  Earth  to  the  radar  site  in  IJK 
coordinates,  while  pUK  is  the  position  vector  from  the  radar  site  to  the  satellite.  The 
components  of  p  can  be  found  in  the  SEZ  frame  quite  easily  using  Figure  4  and  basic 
trigonometry: 
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PsEZ 


—p  cos(cr)  cos(/?) 
p  sin(«)  cos  (/?) 
P  sin(^) 


(48) 


p  can  then  be  converted  to  the  IJK  frame  by  multiplying  it  by  the  inverse  of  the  rotation 
matrix  D: 


P ijk  D  pSEZ 


(49) 


where 


D  = 


sin(Z)  cos(ZST) 
-sin  (LST) 
cos(Z)  cos(ZST) 


sin(Z)sin(ZST)  -cos(Z) 
cos  (LST)  0 

cos(Z)sin(ZST)  sin(Z) 


(50) 


L  represents  the  site’s  geodetic  latitude  and  LST  represents  the  local  sidereal  time.  LST 
is  measured  from  the  vernal  equinox  to  the  radar  site  in  a  counter-clockwise  direction. 
LST  is  also  the  sum  of  the  Greenwich  Sidereal  Time  (GST)  and  the  site’s  longitude. 
Equation  (5 1)  is  used  next  to  calculate  the  position  from  the  center  of  the  Earth  to  the 
radar  site: 


Rs 


xcos(ZST) 

xsin(ZST) 

z 


(51) 


The  quantities  x  and  z  take  into  account  the  fact  that  Earth  is  not  a  perfect  sphere,  but 
rather  an  ellipsoid: 


x  = 


R, 


yj \~e\  sin"  L 


+  H 


cos  L 


(52) 
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^s(i-4) 

yj\  -e^  sin2  L 


+  H 


sinZ 


(53) 


With  Rsite  and  the  values  of  pIJK  for  the  entire  track  in  hand,  the  position  of  the  satellite  in 
IJK  coordinates  is  obtained  for  each  observation.  These  values  are  directly  inserted  into 
the  estimator  as  the  observation  data  (Z.  ]  to  create  an  initial  state  vector  X{t0) . 

Truth  Model 

The  truth  model  is  basically  a  tool  used  to  create  simulated  data  and  an  aid  for 
verifying  the  estimator.  This  model  creates  data  in  a  format  that  a  typical  radar  site 
would  expect:  range,  azimuth,  and  elevation  values.  Due  to  the  specifics  of  this 
research,  the  truth  model  must  be  able  to  simulate  satellite  motion  with  or  without  an 
extra  acceleration  component.  The  truth  model  must  also  simulate  real-world  factors 
such  as  noise.  This  section  explains  how  all  of  these  requirements  were  addressed,  as 
well  as  how  the  truth  model  is  executed  and  validated.  This  section  is  divided  into  three 
sub-sections:  Program  Execution,  Equations,  and  Program  Validation. 

Program  Execution. 

The  truth  model  takes  an  initial  input  for  r,  r,  f,  and  r  in  units  of  km,  km/s, 
km/s",  km/s  respectively.  These  inputs  are  used  in  the  equations  of  motion  to  create  a 
matrix  of  position  and  velocity  vectors  for  a  five-minute  radar  track  for  the  radar  site 
located  at  Ascension,  Atlantic  (another  site  can  be  specified  at  the  beginning  of  the 
program).  Matrices  for  range,  azimuth,  and  elevation  are  created  from  the  position  data 
for  use  in  the  estimator.  The  model  places  the  epoch  time  at  the  middle  of  the  trajectory 
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as  opposed  to  the  beginning.  Output  from  the  truth  model  consists  of  range  (kilometers), 
azimuth  (radians),  elevation  (radians),  time  (in  Julian  days)  for  the  entire  five-minute 
track,  plus  the  site’s  latitude  and  noise  statistics  in  range,  azimuth,  and  elevation. 

The  truth  model  simulation  used  a  start  time  of  13  Sep  2007,  12:00:00  UTC 
(Coordinated  Universal  Time)  with  an  epoch  of  12:02:30  UTC.  Noise  was  added  to  the 
range,  azimuth,  and  elevation  data.  Data  from  the  truth  model  was  output  to  a  file  for 
later  use  in  the  estimator. 

Equations. 

The  first  step  to  simulating  data  and  creating  a  truth  model  was  to  define  the 
dynamics.  Satellite  motion  in  this  case  was  not  modeled  using  two-body  motion,  but  was 
instead  estimated  using  Galileo’s  expertise  on  modeling  projectile  trajectories  with  a 
Taylor  series  expansion.  The  equations  of  motion  below  represent  a  3ld  order  Taylor 
series  approximation: 


r  =  r0  +  v0  A  t  + 1  (A0  +  g0  )A  t2  +  \  (A,  +  So  )At' 

2  6 

(54) 

V  =  V0  +  (A)  +  So  )At  +  -(A)  +  So  )Ar 

(55) 

a  =  (A  +  g0)  +  (20  +  §o)At 

(56) 

a  =  (^o  +  io) 

(57) 

A  discussion  on  why  the  3rd  order  Taylor  series  approximation  was  used  for  the  truth 
model,  and  not  some  other  order,  is  saved  for  the  Analysis  and  Results  chapter.  In  brief, 
the  3ld  order  equations  were  more  accurate  than  the  2nd  order  equations,  and  had  results 
sufficient  enough  to  not  need  the  4th  order  equations. 
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A0  is  the  non-gravitational  acceleration  vector  that  is  to  be  detected.  This 

2 

component  is  assumed  to  be  constant.  The  g  vector  in  this  case  is  not  9.81  m/s  ;  it  is  the 
negative  gradient  of  the  geopotential: 
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geopotential 
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dx 
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dz 
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The  first  term  of  the  geopotential  represents  the  two-body  problem  and  the  second  term 
represents  the  J2  effects  (Wiesel,  2003a:  144).  The  variables  x,  y,  and  z  are  the  Cartesian 


coordinates  of  the  initial  value  of  r  ,  which  can  also  be  written  as  RIJK  .  An  expanded 


view  of  the  partial  derivatives  can  be  seen  in  Equations  (60)-(62). 
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15  z2juJ2Rl 
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(60) 


(61) 


(62) 


The  above  equations  can  be  divided  into  the  g2_body  (Equation  (63))  and  the  gh 
(Equation  (64))  components. 
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Since  the  equations  of  motion  require  the  derivative  of  g ,  it  was  necessary  to 
differentiate.  Keeping  in  mind  that  x,  y,  z  and  r  are  all  functions  of  time,  the  following 
derivatives  for  g2-body  and  the  g  f  were  obtained: 


fjJ 2Rq 


-15x(r®v)  3x  105  xz2(vv)  30  xzz  15  xz 

7  5  q  7  7 

y  y  y  y  y 
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7  s  Q  7 
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3  jUx(r»v)  jUx 
r5  r3 

3  juy(r»v)  juv 

r5  r3 

3  jUz(f»v)  //z 


The  variables  x ,  y ,  and  z  are  the  individual  components  of  the  initial  value  for 
velocity  and  are  written  as  u,  v,  and  w  in  the  MATLAB  code.  The  term  r  •  v  can  be 
written  as  xu+yv+zw.  The  value  of  g  is  the  local  gravity  vector  for  the  given  initial 
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position  and  velocity.  Since  the  values  for  x,  y,  z,  x,y,  and  z  are  taken  from  only  the 
initial  position  and  velocity,  the  value  of  both  g  andg  is  constant  for  the  five-minute 
radar  track. 

After  using  the  equations  of  motion  to  find  the  position  vectors  (in  IJK),  the  data 
was  transfonned  to  range,  azimuth  and  elevation  values  (in  SEZ)  because  this  is  the 
format  of  the  data  that  radar  sites  would  receive.  The  values  for  range,  azimuth,  and 
elevation  change  depending  on  the  location  of  the  radar  site.  The  radar  site  represented  in 
the  truth  model  is  Ascension,  Atlantic.  Table  1  lists  the  specifics  for  this  site. 


Table  1.  Ascension  Atlantic  Radar  Specifics  (Vallado,  2001:242) 


Location 

Noise 

Latitude  (°) 

Longitude  (°) 

Altitude  (m) 

Range (m) 

Azimuth  (°) 

Elevation  (°) 

-7.91 

-14.40 

56.1 

101.7 

0.0248 

0.0283 

The  site  location  was  used  to  detennine  Rsite  (Equations  (5 1  )-(53)).  The  noise 
values  in  Table  1  were  saved  for  use  in  the  estimator.  Figure  5  shows  how  pUK  was 
found  given  Rsite  and  RIJK  .  The  conversion  from  pUK  to  pSE/  was  obtained  by  using 
Equation  (49).  The  values  of  range,  azimuth,  and  elevation  for  each  observation  were 
found  using 


P  =  \IPs+Pe+P- 

(67) 

a  =  tan-1 

(68) 

~Ps 

p  =  suT1  — 

(69) 

P 
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These  radar  values  for  the  entire  time  span  are  saved  in  matrices. 

Radar  sites  will  never  receive  perfect  data;  therefore,  the  truth  model  adds  noise  to 
the  range,  azimuth,  and  elevation  values  to  simulate  real  world  results.  MATLAB’s 
random  number  generator  ( randn )  was  used  to  simulate  noise.  Randn  is  a  pseudo¬ 
random  function  that  creates  noise  with  a  normal  (Gaussian)  distribution.  The  mean  and 
standard  deviation  of  the  noise  can  be  specified  when  using  randn.  Range,  azimuth,  and 
elevation  each  have  a  noise  matrix  associated  with  it.  The  mean  for  these  noise  matrices 
is  zero. 

Program  Validation. 

The  position  and  velocity  matrices  obtained  from  the  equations  of  motion  for  a 
given  initial  position,  velocity,  and  for  both  zero  acceleration  and  its  derivative  were 
compared  to  the  same  initial  inputs  in  a  Satellite  Tool  Kit  (STK)  simulation.  STK  allows 
the  user  to  estimate  an  orbit  using  various  propagation  methods  such  as  two-body,  h,  h 
etc.  Since  the  intent  was  to  verify  how  close  the  MATLAB  model  was  to  reality,  the 
High  Precision  Orbit  Propagator  (HPOP)  was  chosen  to  propagate  the  trajectory.  This 
propagator  takes  into  account  many  elements  such  as:  central  body  gravity,  solar 
radiation  pressure,  third  body  gravity  (sun,  moon)  and  drag,  plus  uses  an  RKF  7(8) 
integrator.  The  initial  conditions  for  the  simulated  orbit  were  written  in  Cartesian 
coordinates  with  position  in  kilometers  and  velocity  in  km/s: 


-6079.6 

1837.9 

-1596.6 


(70) 


30 


(71) 


v 


0 


-2.96 

-5.65 

4.82 


The  magnitude  of  the  position  is  roughly  6549  km.  Taking  into  account  that  the 
radius  of  the  Earth  is  approximately  6738  km,  the  simulated  orbit  has  an  altitude  of  only 
171  km.  This  altitude  is  a  very  low  LEO  and  most  objects  do  not  last  long  at  such  a  low 
altitude  due  to  atmospheric  drag.  However,  this  altitude  could  be  quite  realistic  for  a 
tethered  satellite  or  ballistic  missile. 

The  results  for  position  and  velocity  between  STK  and  MATLAB  were  very 
close,  but  changed  in  accuracy  depending  on  which  order  of  the  Taylor  series 
approximation  was  used.  A  visual  comparison  of  the  accuracy  of  the  2nd,  3rd,  and  4th 
order  approximations  can  be  seen  in  Chapter  IV.  Table  2  shows  the  root  mean  square 
(RMS)  error  for  the  x,  y,  z,  x ,  y ,  and  z  components  of  position  and  velocity. 


Table  2.  STK  vs.  MATLAB  RMS  Error  (3rd  Order  EOM) 


Position  (km) 

Velocity  (km/s) 

X 

y 

z 

X 

y 

z 

0.096835 

0.034235 

0.029878 

0.002951 

0.00099 

0.000866 

The  values  in  Table  2  were  calculated  using  Equation  (72): 


RMS  = 


z 


(72) 
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The  component  et  is  the  individual  magnitudes  of  error  between  the  MATLAB  and  STK 

values  and  N  is  the  number  of  data  points.  In  this  case,  N=300,  one  for  each  second 
during  the  five-minute  track. 

STK  was  also  used  to  validate  the  local  gravity  vector  by  finding  the  satellite’s 
acceleration  component  at  epoch  in  the  STK  simulation.  Given  the  above  initial 
conditions,  the  local  gravity  vector  calculated  by  the  truth  model  is 


g  = 


0.00863715 

-0.0026111 

0.00227524 


km/s2 


(73) 


This  is  extremely  close  to  STK’s  acceleration  value  at  epoch: 


g  = 


0.008638 

-0.002610 

0.002275 


km/s2 


(74) 


The  range,  azimuth,  and  elevation  portion  of  the  code  was  validated  with  the  STK 
simulation  as  well  as  with  code  from  the  previous  section. 


Estimator 

Execution. 

The  estimator  reads  in  an  input  file  that  contains:  elevation,  range,  azimuth,  time 
(in  Julian  days),  site  latitude,  and  noise  statistics  for  range,  azimuth,  and  elevation.  The 
truth  model  used  units  of  kilometers,  radians  and  Julian  days  for  the  various  components 
of  output  data;  however,  radar  sites  typically  get  their  data  in  units  of  meters,  degrees, 
and  time  (in  year,  day  number,  hour,  minutes,  and  seconds).  Therefore,  the  estimator  has 
an  option  where  the  radar  data  can  be  converted  into  the  kilometers,  radians,  and  Julian 
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day  format.  The  data  is  then  converted  to  position  vectors  using  the  process  outlined  in 
Section  3.2.  The  estimator  uses  linear  least  squares  to  determine  the  estimate  of  the  state 
vector  at  epoch  to  include:  r0 ,  v0 ,  a0 ,  and  a0 .  The  estimator  also  produces  the 
covariance  of  the  state  vector.  Since  this  is  a  linear  system,  there  is  no  need  to  compute 
residuals  and  iterate.  There  is  also  no  need  for  an  initial  guess  of  the  state  vector. 

Equations. 

The  state  vector  for  this  problem  was  defined  as 

r 

v 


X  = 


(75) 


a 

a 


The  initial  state  vector  looks  like 


7) 

Vo 

(4)  +  g0) 

a0 

_U)+£o)J 

A_ 

(76) 


With  this  definition  of  the  state  vector  and  the  dynamics  of  the  system  already  known,  a 
closed  form  solution  of  the  state  transition  matrix  was  obtained  using  Equation  (35)  from 
Chapter  II: 


33 


dr 

dr 

dr 

dr 

r)v 

Ut  0 

dv0 

da0 

da0 

dv 

dv 

dv 

dv 

ur  0 

dv0 

da0 

da0 

da 

da 

da 

da 

df0 

1 l>0 

1  rD 

1 

1  rD 

dd0 

da 

da 

da 

da 

r)v 
ur  0 

dv0 

da0 

da0 

(77) 


Each  component  is  a  3x3  matrix,  which  makes  the  state  transition  matrix  a  square 
12x12  matrix.  The  diagonal  terms  are  the  identity  matrix.  All  terms  to  the  left  of  the 

0—  3v 

diagonal  are  the  null  matrix.  The  terms  — and  — ^  are  the  identity  matrix 

dv0  da0  da0 


dr  dv  t~  dr 

multiplied  by  time.  The  terms  — and  — ^  are  the  identity  times  — ,  while  — ^  is  the 


da, 


da 


da,, 


identity  times  — .  A  multi-dimensional  array  for  the  state  transition  matrix,  with  a  one- 
6 


second  time  step,  within  the  five-minute  track  is  created. 

In  the  quest  to  simplify  the  observation  relation  G  and  to  linearize  the  data,  the 
observed  data  vector  consisting  of  range,  azimuth,  and  elevation  components  is  converted 
to  pseudo-data  (Wiesel,  2003b:94-95): 


Z’  =  f 

Z'  =  G(X,t)  =  {I,<p,<p,<p)X 


(78) 


In  this  equation,  /  is  the  identity  matrix  and  (p  is  the  null  matrix.  This  creates  a 
simplified  G  function  with  its  linearization  H  also  simplified  and  not  a  function  of  time: 


34 


"1  0000000000  0" 

~\fy 

H  =  ~=  010000000000  (79) 

dx 

001000000000 

After  H  is  obtained,  the  observation  matrix  (T)  can  be  obtained  (see  Equation  (39)). 

The  use  of  pseudo-data  instead  of  range,  azimuth,  and  elevation  creates  a  Q 
matrix  that  is  no  longer  constant  (Wiesel,  2003b:95).  Therefore,  the  values  of  Q  must  be 
calculated  for  every  position.  Wiesel  shows  that  these  values  are  easily  obtained  through 
a  simple  rotation: 

Q'  =  JQJT  (80) 

with  the  original  covariance  written  as 

Perror  0  0 

Q=  0  a;rmr  0  (81) 

_  o  o  ffmr_ 

The  values  of  perror ,  aerror ,  and  f3error  are  obtained  from  Table  1. 

The  Jacobian  is  obtained  by 

J  =  DlK  (82) 

where 

-  cos  /?  cos  a  /?  cos/?  sin  or  /?  sin/?  cos  or 
K=  cos/?  sin  a  p  cos  ftcosa  -p  sin  /?  sin  a  (83) 

sin  /?  0  p  cos  /? 

and  D1  is  the  inverse  of  Equation  (50). 

All  of  the  necessary  matrices  needed  to  find  the  state  vector  at  epoch  and  its 
covariance  have  been  found.  Since  least  squares  is  a  batch  process,  the  least  squares 
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equations  obtained  in  Chapter  II,  Equations  (27)  and  (28),  can  be  reformatted  and  written 


as 


f  N 


s-1 


Zrre;‘7; 


VH 


(84) 


xih)  =  Pxfr?Q;%  (85) 

i=i 

in  order  to  save  computer  space  and  ensure  a  quicker  processing  time.  In  this  case,  the 
product  of  the  covariance  and  state  vector  are  summed  for  a  five-minute  radar  track  from 
t=  -150  to  t=T49  seconds.  After  the  state  vector  at  epoch  is  obtained,  any  state  vector 
thereafter  is  obtained  simply  by  using  Equation  (36). 

Validation. 

The  truth  model  was  used  to  validate  the  estimator.  Various  initial  positions  and 
velocities  without  non-gravitational  acceleration  were  run  through  the  truth  model  and 
the  estimator  to  verify  that  the  estimated  state  was  within  la  of  the  true  state.  The 
estimated  state  and  its  covariance  were  also  validated  using  the  Monte  Carlo  method. 

The  Monte  Carlo  method  was  used  to  ensure  that  the  covariance  was  representative  of  the 
uncertainty  in  the  data.  An  in-depth  look  at  the  Monte  Carlo  method  and  its  results  is 
included  in  Chapter  IV. 


Summary 

Numerous  MATLAB  routines  were  created  in  order  to  develop  the  truth  model 
(Appendix  B)  and  least  squares  estimator  (Appendix  C).  Given  an  initial  position  and 
velocity,  the  truth  model  develops  a  five-minute  track  of  data  and  converts  it  to  the 
familiar  range,  azimuth,  and  elevation  format.  The  least  squares  estimator  takes  both  the 
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radar  data  and  instrument  covariance  in  order  to  estimate  the  epoch  value  of  the  state 
vector  (to  include  r,  r,  f,  and  r  )  and  the  covariance  of  the  state.  A  flowchart  that  depicts 
this  process,  from  the  initial  input  data  to  the  final  estimate  of  the  state  and  its  covariance, 
is  located  in  Appendix  D.  The  magnitude  of  the  additional  acceleration  and  its 
covariance  were  the  components  of  interest  in  this  research.  A  more  thorough  analysis  of 
the  acceleration  and  its  covariance  is  discussed  in  the  next  chapter. 
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IV.  Analysis  and  Results 


Chapter  Overview 

The  purpose  of  this  chapter  is  to  analyze  the  various  sets  of  equations  of  motion 
used  in  the  truth  model  to  model  a  space  object  in  LEO;  and  to  detennine  if  one  of  these 
sets  of  equations  of  motion  combined  with  a  linear  least  squares  estimator  can 
satisfactorily  detect  non-gravitational  acceleration  with  statistical  accuracy.  Section  4.2 
explores  the  accuracy  of  the  various  sets  of  equations  of  motion  compared  to  an  STK 
simulation  of  the  same  orbit.  Sections  4.3  and  4.4  present  several  test  cases  with  various 
initial  conditions  used  to  analyze  how  well  the  truth  model  and  estimator  function.  These 
test  cases  help  to  detennine  the  degree  of  non-gravitational  acceleration  that  can  be 
adequately  detected.  Section  4.5  addresses  the  investigative  questions  that  were  posed  in 
Chapter  I.  Section  4.6  summarizes  the  main  discoveries  of  the  research. 

STK  Simulation  vs.  MATLAB  Model 

The  initial  conditions  stated  in  Equations  (70)  and  (71)  were  used  to  create  a  STK 
simulation  for  use  as  a  baseline  model.  The  2nd,  3ld,  and  4th  order  Taylor  series 
approximations  were  compared  to  the  STK  simulation  to  detennine  which  equations 
modeled  a  space  object  accurately  enough.  The  same  initial  conditions  were  input  into 
the  truth  model.  For  this  analysis,  both  the  non-gravitational  acceleration  and  its 
derivative  were  assumed  to  be  zero.  After  obtaining  the  position  and  velocity  values 
from  the  truth  model  for  the  five-minute  radar  track,  these  values  were  compared  to  the 
STK  values.  The  amount  of  enor  between  the  STK  values  and  MATLAB  values  was 
calculated  using 
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error  =  observed  -  expected  =>  MATLAB  -  STK 


(86) 


The  values  of  error  for  both  position  and  velocity  were  graphed  using  Excel  in 
order  to  get  a  visual  idea  of  the  accuracy  of  the  equations.  The  following  sub-sections  go 
into  detail  about  the  accuracy  of  the  various  orders  of  the  equations  of  motion. 

2nd  Order. 

The  2nd  order  equations  of  motion  were  the  initial  equations  tested  for  the  truth 
model.  Figure  6  shows  the  relative  error  in  position  compared  to  the  STK  simulation.  As 
would  be  expected,  the  error  is  minimal  around  epoch  but  grows  as  time  moves  on.  The 
amount  of  error  at  these  other  times  is  not  particularly  ideal. 


Figure  6.  STK  vs.  MATFAB  Position  Error  (2nd  Order) 


Figure  7  shows  the  error  in  velocity  compared  to  the  STK  simulation.  As  seen  in 
Figure  6,  the  error  is  zero  at  epoch  but  grows  significantly  as  time  moves  away.  Fooking 
at  the  z  velocity  line,  the  greatest  magnitude  in  error  is  at  t=0  seconds  where  velocity  is 


39 


roughly  0.08  km/s  or  80  m/s.  The  error  magnitude,  0.53  m/s1 ,  was  obtained  by  dividing 
the  velocity  by  half  the  observation  time 

80  m/s 


150^ 


=  0.53  m/s 


(87) 


This  value  is  potentially  the  amount  of  acceleration  that  could  go  unnoticed  due  to  the 
level  of  error  in  the  equations  of  motion.  The  magnitude  of  error  in  this  model  is  too  high 
for  the  2nd  order  equations  of  motion  to  be  of  any  use.  The  model  must  be  able  to  detect 
accelerations  in  the  cm/s  or  possibly  even  pm/s“  range;  therefore,  the  equations  of 
motion  must  be  expanded  out  to  obtain  better  accuracy. 


- x  velocity 

- y  velocity 

z  velocity 


Figure  7.  STK  vs.  MATLAB  Velocity  Error  (2nd  Order) 

3rd  Order. 

The  equations  of  motion  were  expanded  out  to  the  3ld  order  to  obtain  a  higher 
degree  of  accuracy.  Figure  8  shows  the  amount  of  position  error  for  the  3ld  order 
equations  of  motion.  As  was  seen  in  the  previous  graphs,  the  amount  of  error  at  epoch  is 
minimal  and  then  grows  as  time  increases.  Unlike  the  previous  graphs,  however,  the 
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magnitude  of  error  is  much  less.  The  greatest  magnitude  of  position  error  in  Figure  6  is 
roughly  4.5  km;  whereas,  the  greatest  error  in  Figure  8  is  only  0.3  km.  This  is  a  15  times 
improvement. 


- x  position 

- y  position 

z  position 


Figure  8.  STK  vs.  MATLAB  Position  Error  (3rd  Order) 


So  far,  the  use  of  the  3ld  order  equations  of  motion  seems  to  be  yielding  better  results. 
Figure  9  shows  the  velocity  error  for  the  3ld  order  equations  of  motion.  These  results  are 
also  much  better  than  that  of  the  2nd  order.  The  greatest  magnitude  of  error  at  t=0  is 
0.008  km/s.  Dividing  this  value  by  150  s  yields 

8  m/s 


150.9 


=  0.053  m/s 


(88) 


Therefore,  the  amount  of  undetected  acceleration  that  could  be  present  in  the  3rd  order 


equations  of  motion  is  0.053  m/sz.  This  result  is  10  times  better  than  the  2nd  order 


equations  of  motion.  If  the  3ld  order  yielded  much  better  results,  it  seems  only  reasonable 


that  expanding  to  the  4th  order  would  obtain  an  even  higher  level  of  detection. 
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Figure  9.  STK  vs.  MATLAB  Velocity  Error  (3rd  Order) 


4th  Order. 

The  4th  order  equations  of  motion  yielded  graphs  that  were  a  bit  different  than 
those  seen  in  the  previous  figures.  Figure  10  represents  the  position  error  for  the  4th  order 
expansion.  Just  like  the  previous  graphs,  there  is  zero  error  at  epoch;  but  instead  of 
growing  exponentially  thereafter,  the  graph  curves  again  at  t=250  seconds.  The  amount 
of  error  present  is  also  significantly  less  than  that  of  the  previous  graphs. 
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Figure  10.  STK  vs.  MATLAB  Position  Error  (4th  Order) 


The  velocity  graph  seen  in  Figure  1 1  also  yielded  much  better  results.  The 

2 

maximum  error  at  t=0  seconds  is  roughly  0.0004  km/s  .  Following  the  same  process  seen 
in  Equations  (87)  or  (88),  the  amount  of  undetected  acceleration  is  0.00267  m/s  .  The 
level  of  detection  is  roughly  20  times  better  than  what  was  seen  using  the  3ld  order 
equations. 
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Figure  1 1 .  STK  vs.  MATLAB  Velocity  Error  (4th  Order) 


Besides  the  visual  comparison  of  the  accuracy  of  the  different  Taylor  series 
expansions,  there  is  the  numerical  approach.  The  RMS  error  for  all  components  in 
position  and  velocity  was  calculated  using  Equation  (72)  for  the  2nd,  3rd,  and  4th  order 
equations  of  motion.  Table  3  lists  the  results.  As  is  expected,  the  RMS  error  decreases  as 
the  order  used  increases.  The  4th  order  may  yield  the  best  results;  however,  the  3ld  order 
results  are  also  quite  viable  compared  to  the  fairly  inadequate  results  seen  from  the  2nd 
order. 


Table  3.  STK  vs.  MATLAB  RMS  Error 


Position  (km) 

Velocity  (km/s) 

X 

y 

z 

X 

y 

z 

2nd  Order 

0.91264 

1.69929 

1.45603 

0.02163 

0.04016 

0.03436 

3  rd  Order 

0.096835 

0.034235 

0.029878 

0.002951 

0.00099 

0.000866 

4th  Order 

0.003297 

0.004272 

0.003524 

8.88E-05 

0.000121 

0.000102 
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The  question  is,  “Which  order  of  equations  of  motion  is  accurate  enough  for  the 
purposes  of  this  research?”  The  position  and  velocity  data  obtained  from  the  truth  model 
for  the  3rd  order  equations  of  motion  was  processed  through  the  least  squares  estimator  to 
obtain  the  state  vector  and  its  covariance.  The  4th  order  data  was  also  processed  through 
the  estimator.  The  4th  order  equations  of  motion  produced  a  15x1  state  vector  while  the 
3ld  order  equations  of  motion  produced  a  12x1  state  vector,  as  seen  in  Table  4. 


Tab 


e  4.  MATLAB  Estimated  State  Vector  with  No  N 


Variable 

4th  Order 

3rd  Order 

x  (km) 

-6079.6 

-6079.6 

y  (km) 

1837.9 

1837.9 

z  (km) 

-1596.6 

-1596.6 

X  (km/s) 

-2.96 

-2.96 

y  (km/s) 

-5.65 

-5.65 

z  (km/s) 

4.82 

4.82 

ax  (km/s2) 

0.008637157 

0.008637157 

av  (km/s2) 

-0.002611065 

-0.002611065 

az  (km/s2) 

0.002275235 

0.002275235 

ax  (km/s3) 

4.2799E-06 

4.2799E-06 

av  (km/s3) 

8.00425E-06 

8.00425E-06 

b_  (km/s3) 

-6.84906E-06 

-6.84906E-06 

tix  (km/s4) 

-1.41524E-08 

tiY  (km/s4) 

4.41528E-09 

d_(km/s4) 

-3.90375E-09 

oise 


The  estimates  of  the  state  vector  for  both  the  3rd  and  4th  order  equations  of  motion 
in  a  noiseless  scenario  were  nearly  identical.  Differences  were  seen,  however,  in  their 
covariance.  For  a  given  covariance  matrix,  the  top  left  diagonal  tenn  represents  the 
variance  of  the  x  position  tenn.  The  second  tenn  is  the  variance  of  the  y  position  and  so 
forth  all  the  way  down  to  the  bottom  right  diagonal  tenn,  which  in  the  3rd  order  case,  is 
the  last  component  of  a  .  Table  5  lists  the  variances  obtained  for  the  noiseless  state. 
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Table  5.  MATLAB  Variance  for  Noiseless  Scenario 


Variable 

4th  Order 

3rd  Order 

x  (km) 

0.002266 

0.001569 

y  (km) 

0.001178 

0.000944 

z  (km) 

0.000549 

0.000482 

X  (km/s) 

7.1E-07 

6.7E-07 

y  (km/s) 

4.02E-07 

2.55E-07 

Z  (km/s) 

2.32E-07 

1.27E-07 

ax  (km/s2) 

8.98E-10 

8.99E-11 

av  (km/s2) 

3.27E-10 

4.89E-11 

a:  (km/s2) 

1.64E-10 

3.18E-11 

ax  (km/s3) 

1.38E-13 

1.3E-13 

av  (km/s3) 

7.51E-14 

4.63E-14 

a _  (km/s3) 

4.9E-14 

2.77E-14 

dx  (km/s4) 

3.63E-16 

av  (km/s4) 

1.29E-16 

a  ,  (km/s4) 

7.7E-17 

The  results  from  Tables  4  and  5  use  data  from  a  noiseless  environment,  which, 
unfortunately,  does  not  accurately  portray  reality.  Gaussian  noise  was  added  to  the 
position  and  velocity  data  for  both  the  3rd  and  4th  order  simulations  to  obtain  new 
estimates  of  the  state  and  covariance.  Table  6  lists  the  estimates  of  the  states.  As 
expected,  these  new  values  are  slightly  different  than  the  values  with  no  noise. 
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Table  6.  MATLAB  Estimated  State  Vector  with  Noise 


Variable 

4th  Order 

3rd  Order 

x  (km) 

-6079.585554 

-6079.5851 

y  (km) 

1837.895358 

1837.8787 

z  (km) 

-1596.596015 

-1596.5725 

X  (km/s) 

-2.959402644 

-2.9605417 

y  (km/s) 

-5.650467507 

-5.6498852 

z  (km/s) 

4.819993529 

4.8202066 

ax  (km/s2) 

0.008611162 

0.0086256 

ay  (km/s2) 

-0.002594883 

-0.0026044 

a.  (km/s2) 

0.00226863 

0.0022762 

Clx  (km/s3) 

4.07297E-06 

4.513E-06 

av  (km/s3) 

8.18258E-06 

7.943E-06 

d_(km/s3) 

-6.85554E-06 

-6.763E-06 

ax  (km/s4) 

2.76991E-09 

a  v  (km/s4) 

-6.81815E-09 

a.  (km/s4) 

1.0762E-09 

As  with  any  estimator,  it  is  necessary  to  ensure  that  the  results  are  unbiased  and 
that  the  covariance  matrix  accurately  reflects  the  amount  of  uncertainty  in  the  data.  Thus 
far,  the  results  of  the  covariance  matrix  have  not  been  validated.  A  technique  called 
Monte  Carlo  analysis  is  often  used  to  validate  the  function  of  the  estimator. 

Monte  Carlo  Analysis. 

There  are  several  steps  required  in  order  to  use  the  Monte  Carlo  method.  For  a 
given  trajectory,  the  truth  model  and  estimator  must  produce  N  number  of  data  sets. 

Each  data  set  has  different  noise,  but  with  the  same  mean  and  standard  deviation. 
Ultimately,  this  produces  slightly  different  estimates  of  the  state  vector.  These  N 
estimates  are  used  to  confirm  that  the  estimator  is,  “  i)  on  the  average  unbiased,  ii)  that 
the  average  estimate  is  the  true  value,  and  iii)  that  the  output  covariance  is  actually 
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representative  of  the  uncertainty  in  the  estimate”  (Wiesel,  2003b:  138).  As  long  as  the 
estimator  is  unbiased,  the  true  state  should  be  obtainable  using  Equation  (89): 

1  N  _ 

Unlike  in  reality,  the  true  state  ( A0 )  is  known  because  it  was  chosen  in  the  truth  model. 


The  variable  Xt  is  the  different  estimates  of  the  state.  Using  this  same  method,  the 
covariance  of  the  state  should  be 


/’“tfIU-UU-U"  (*>> 

N  ;= l 

Up  to  a  certain  point,  increasing  the  number  of  data  sets  yields  more  accurate  results. 

Ten  data  sets  were  obtained  for  the  test  case.  Table  7  displays  data  sets  1-5  where  the 
results  are  listed  in  the  order  XT  =  |V  v  a  a  a~^. 


Table  7.  State  Vector  Data  Sets  1-5  (4th  Order) 


Set  1 

Set  2 

Set  3 

Set  4 

Set  5 

x  (km) 

-6079.607723 

-6079.669544 

-6079.627863 

-6079.585554 

-6079.588185 

y  (km) 

1837.866979 

1837.957821 

1837.940816 

1837.895358 

1837.894576 

z  (km) 

-1596.559154 

-1596.628007 

-1596.627217 

-1596.596015 

-1596.598784 

X  (km/s) 

-2.960495055 

-2.957693797 

-2.960761762 

-2.959402644 

-2.959472693 

y  (km/s) 

-5.650189915 

-5.650678485 

-5.648872485 

-5.650467507 

-5.650354919 

z  (km/s) 

4.820082613 

4.819354213 

4.81885531 

4.819993529 

4.820127364 

ax  (km/s2) 

0.008648758 

0.008666269 

0.008630577 

0.008611162 

0.008616777 

av  (km/s2) 

-0.00259395 

-0.002621308 

-0.00263258 

-0.002594883 

-0.002595558 

a,  (km/s2) 

0.002266441 

0.002273629 

0.002259587 

0.00226863 

0.002270821 

ax  (km/s3) 

4.47891E-06 

3.49766E-06 

4.72945E-06 

4.07297E-06 

3.92502E-06 

d  (km/s3) 

8.08021E-06 

8.33767E-06 

7.45924E-06 

8.18258E-06 

7.99297E-06 

b,  (km/s3) 

-6.70901E-06 

-6.744 18E-06 

-6.40824E-06 

-6.85554E-06 

-6.90807E-06 

tix  (km/s4 

-2.94E-08 

-3.25441E-08 

-1.77783E-09 

2.76991E-09 

4.34767E-09 

bv  (km/s4 

-2.08876E-09 

5.79647E-09 

1.59187E-08 

-6.81815E-09 

-1.81219E-09 

b_  (km/s4) 

1.59348E-09 

-4.27323E-09 

1.14212E-08 

1.0762E-09 

-4.33308E-10 
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Plugging  the  ten  data  sets  into  Equations  (89)  and  (90)  yields  an  estimate  of  the 
state  vector  and  its  covariance.  The  diagonal  terms  of  the  Monte  Carlo  estimate  of  the 
covariance  matrix  for  the  4th  order  are  listed  below  in  Table  8.  It  is  important  to 
remember  from  Chapter  II  that  the  square  root  of  the  variance  determines  the  standard 
deviation.  Of  notable  interest  are  the  results  in  columns  4  and  5.  Column  4  displays  the 
Monte  Carlo  estimation  of  both  the  true  state  and  its  standard  deviation.  Careful  analysis 
shows  that  the  true  value  of  the  state  is  within  the  bounds  of  the  estimated  value  except  in 
the  a  components.  The  estimated  values  and  the  true  values  are  quite  off.  Also,  the 
magnitude  of  the  estimated  a  component  is  on  order  of  10'1 1  km/s3  with  a  much  larger 

standard  deviation  on  order  of  10'  km/s  . 

Essentially,  the  estimator  is  incapable  of  properly  estimating  these  small 
magnitudes.  Evidence  of  this  can  be  verified  above  in  Table  7.  All  five  data  sets  have 
vastly  different  values  for  the  a  components.  Since  the  estimator  is  unable  to  accurately 
estimate  the  a  components,  it  is  inefficient  to  use  the  higher  order  equations.  It  seems 
quite  reasonable  to  go  down  a  level  of  accuracy,  and  simply  use  the  3ld  order  equations. 
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Table  8.  Monte  Carlo  Estimated  State  Vector  Values  (4th  Order) 


Variable 

Variance 

Sigma 

Estimated  Value 

True  Value 

x  (km) 

0.0023495 

0.048471641 

-6079.606304  ±0.048471641 

-6079.6 

y  (km) 

0.002409 

0.049081565 

1837.916153  ±0.049081565 

1837.9 

z  (km) 

0.00078343 

0.02798982 

-1596.608641  ±0.02798982 

-1596.6 

X  (km/s) 

1.2272E-06 

0.001107791 

-2.959794471  ±0.001107791 

-2.96 

y  (km/s) 

5.4565E-07 

0.000738681 

-5.650153622  ±  0.000738681 

-5.65 

z  (km/s) 

2.1741E-07 

0.000466272 

4.819836215  ±0.000466272 

4.82 

ax  (km/s2) 

1.2268E-09 

3.50257E-05 

0.00863 1874  ±  3.50257E-05 

0.008637157 

av  (km/s2) 

5.6799E-10 

2.38325E-05 

-0.002609078  ±2.38325E-05 

-0.002611065 

O,  (km/s2) 

2.1518E-10 

1.4669E-05 

0.0022729  ±  1.4669E-05 

0.002275235 

ax  (km/s3) 

2.0364E-13 

4.51265E-07 

4.21247E-06  ±  4.51265E-07 

4.2799E-06 

av  (km/s3) 

1.0374E-13 

3.22087E-07 

8.0802E-06±  3.22087E-07 

8.00425E-06 

d,  (km/s3) 

4.3888E-14 

2.09495E-07 

-6.8051E-06±  2.09495E-07 

-6.84906E-06 

tix  (km/s4) 

6.9845E-16 

2.64282E-08 

-7.10425E-09  ±  2.64282E-08 

-1.41524E-08 

dv  (km/s4) 

1.5985E-16 

1.26432E-08 

2.8870  IE- 1 1  ±  1.26432E-08 

4.41528E-09 

a,  (km/s4) 

8.2957E-17 

9.10807E-09 

-1.37777E-09  ±  9.10807E-09 

-3.90375E-09 

The  process  used  for  obtaining  the  4th  order  data  sets  was  also  used  to  obtain  the 
3ld  order  data  sets.  Table  9  displays  data  sets  1-5.  By  inspection,  it  is  apparent  that  the 
estimates  in  each  row  all  have  values  that  are  quite  close. 


Table  9.  State  Vector  Data  Sets  1-5  (3rd  Order) 


Set  1 

Set  2 

Set  3 

Set  4 

Set  5 

x  (km) 

-6079.585097 

-6079.644248 

-6079.637398 

-6079.621801 

-6079.619295 

y  (km) 

1837.878714 

1837.955713 

1837.93543 

1837.901333 

1837.910241 

z  (km) 

-1596.572505 

-1596.631851 

-1596.62999 

-1596.583491 

-1596.598289 

X  (km/s) 

-2.960541698 

-2.957667187 

-2.96116703 

-2.957978792 

-2.95971252 

y  (km/s) 

-5.649885182 

-5.650481368 

-5.649712807 

-5.650918108 

-5.650543327 

z  (km/s) 

4.820206618 

4.819342736 

4.819486692 

4.820305664 

4.820406233 

ax  (km/s2) 

0.008625584 

0.008638831 

0.008647061 

0.008642268 

0.008646995 

a  (km/s2) 

-0.002604419 

-0.002619199 

-0.002618182 

-0.002611856 

-0.002608052 

fl,  (km/s2) 

0.00227618 

0.002274249 

0.00227742 

0.002271089 

0.002280282 

ax  (km/s3) 

4.51282E-06 

3.49377E-06 

4.90508E-06 

3.60324E-06 

4.07412E-06 

a  (km/s3) 

7.94287E-06 

8.2494E-06 

7.83093E-06 

8.33482E-06 

8.0473E-06 

d,  (km/s3) 

-6.76297E-06 

-6.73871E-06 

-6.69302E-06 

-7.042E-06 

-6.88487E-06 
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Table  10  displays  the  results  for  the  estimated  state  and  its  covariance  using  the 
3ld  order  equations  of  motion  data  sets.  For  every  component  of  the  state  vector,  the  true 
state  is  within  the  bounds  of  the  estimated  value  and  standard  deviation. 


Table  10.  Monte  Carlo  Estimated 

l  State  Vector  Values  (3rd  Order) 

Variable 

Variance 

Sigma 

Estimated  Value 

True  Value 

x  (km) 

0.00078788 

0.0280692 

-6079.617422  ±  0.0280692 

-6079.6 

y  (km) 

0.0016537 

0.04066571 

1837.92279  ±0.04066571 

1837.9 

z  (km) 

0.0006166 

0.024831432 

-1596.609597±  0.024831432 

-1596.6 

X  (km/s) 

1.69990000E-06 

0.001303802 

-2.959698749  ±  0.001303802 

-2.96 

y  (km/s) 

3.12090000E-07 

0.00055865 

-5.650207329  ±0.00055865 

-5.65 

Z  (km/s) 

1.25700000E-07 

0.000354542 

4.819932336  ±0.000354542 

4.82 

ax  (km/s2) 

1.02870000E-10 

1.01425E-05 

0.008643123  ±  1.01425E-05 

0.008637157 

aY  (km/s2) 

9.94220000E-1 1 

9.97105811E-06 

-0.002616  ±  9.9710581 184E-06 

-0.002611065 

(km/s2) 

4.20220000E-1 1 

6.48243781E-06 

0.00227658  ±  6.482437813E-06 

0.002275235 

ax  (km/s3) 

2.58550000E-13 

5.08478121E-07 

4.19095E-06  ±  5.08478121E-07 

4.2799E-06 

av  (km/s3) 

5.42150000E-14 

2.32841 147E-07 

8.09606E-06  ±  2.3284 1 147E-07 

8.00425E-06 

d,  (km/s3) 

2.1 1380000E-14 

1.45389133E-07 

-6.83767E-06±  1.45389133E-07 

-6.84906E-06 

The  variances  obtained  in  Tables  8  and  10  have  different  magnitudes  than  their 
respective  variances  in  Table  5.  Fortunately,  there  is  a  plausible  explanation  for  this. 

The  Monte  Carlo  analysis  is  a  weighted  average,  not  an  exact  answer.  The  uncertainty  in 

a  weighted  average  drops  off  proportional  to  — J=  .  In  this  case,  since  there  are  only  ten 

y/N 

data  sets,  the  uncertainty  in  the  Monte  Carlo  covariance  matrix  is  roughly  32%. 
Producing  more  data  sets  would  decrease  the  uncertainty  but  there  will  come  a  point 
where  an  enormous  amount  of  data  sets  is  required  to  improve  results  by  only  a  fraction 
of  a  percentage. 

Based  on  the  above  results,  the  3rd  order  equations  of  motion  appear  quite 
efficient  at  modeling  an  object  in  orbit  over  very  short  arcs.  The  above  results  also 
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support  the  conclusion  that  the  estimator  is  sufficient  at  obtaining  the  state  and  its 
covariance  for  the  3ld  order  equations  of  motion.  The  following  sections  will  examine 
how  well  the  3ld  order  equations  of  motion  and  the  computer  model  are  able  to  estimate 
orbits  with  different  initial  conditions  for  position  and  velocity,  as  well  as  estimate  orbits 
with  various  magnitudes  of  non-gravitational  acceleration  present. 

MATLAB  Simulations  with  Zero  Non-Gravitational  Acceleration 

Equations  (70)  and  (71)  contain  the  initial  conditions  of  the  orbit  that  was  used  to 
validate  the  accuracy  of  the  truth  model  to  that  of  STK.  The  estimator  also  proved 
capable  of  estimating  these  initial  conditions  given  a  five-minute  radar  track  of  data.  It  is 
important,  however,  to  ensure  that  the  estimator  is  capable  of  estimating  the  state  given 
different  values  for  the  initial  position  and  velocity.  A  few  cases  with  different  initial 
conditions  were  tested.  The  first  case  was  comprised  of  the  following  components 

'  1611  ' 

r0  =  -1756  km  (91) 

6100 

7 

v0=  0.4  km/s  (92) 

3.84 

The  second  case  consisted  of 

"  6250  " 

r0  =  500  km  (93) 

-1891 
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km/s 


(94) 


Both  of  these  cases  are  in  LEO  with  roughly  a  170  km  altitude  and  a  velocity  magnitude 
of  about  7.99  km/s.  Gaussian  noise  was  added  to  the  data  and  the  following  estimates  in 
Table  1 1  were  obtained  for  a  single  run. 


Table  11.  Estimated  State  Vectors 


Case  A 

Case  B 

x  (km) 

1611.22 

6250.073 

y  (km) 

-1755.86 

500.7342 

z  (km) 

6099.83 

-1891.45 

X  (km/s) 

6.99748 

1.700548 

y  (km/s) 

0.40367 

-5.00128 

Z  (km/s) 

3.844742 

6.003245 

ax  (km/s2) 

-0.00234 

-0.00889 

ay  (km/s2) 

0.00249 

-0.00079 

a  „  (km/s2) 

-0.00855 

0.002731 

ax  (km/s3) 

-3.2E-06 

-4.6E-06 

ay  (km/s3) 

-7.2E-06 

6.42E-06 

a.  (km/s3) 

1.33E-05 

-7.3E-06 

Based  on  the  above  table,  it  appears  that  the  estimator  is  quite  capable  of  detennining  the 
state  with  different  initial  position  and  velocity  values. 


MATLAB  Simulations  with  Non-Gravitational  Acceleration  Present 

Given  the  equations  of  motion  listed  in  Equations  (54)-(57),  the  truth  model  and 
estimator  are  limited  in  the  magnitude  of  non-gravitational  acceleration  that  can  be 
detected.  In  order  to  detennine  this  magnitude,  it  is  necessary  to  solve  the  linear  first- 
order  differential  equation: 
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(95) 


r  =  —A0t2  +  —  A0t3 

2  o  6  o 

Assuming  that  the  radar  sensor  has  a  range  error  of  100  m,  the  non-gravitational 
acceleration  acting  on  the  space  object  must  move  the  object  more  than  100  m  between 
its  first  and  last  point  in  order  to  statistically  validate  that  acceleration  is  present; 
therefore,  the  vector  r  is  100  m.  Solving  for  A0  yields  the  general  solution 


-  600  C 

An  =^^  +  -r 


(96) 


Time  is  the  length  of  the  observation,  which  in  this  case  is  five  minutes.  A  particular 
solution  for  A0  could  be  obtained  if  some  initial  conditions  were  known.  In  reality, 

however,  the  exact  value  of  the  constant  (C)  is  unobtainable.  If  C=0,  the  magnitude  of 
non-gravitational  acceleration  required  for  detection  given  the  equations  of  motion  is 
0.667  cm/s".  This  value  is  the  smallest  allowable  magnitude  given  the  aforementioned 
conditions.  Given  this  requirement,  all  test  cases  used  in  this  research  have  a  magnitude 
greater  than  or  equal  to  0.667  cm/s". 

Monte  Carlo  Approach. 

Six  test  cases  were  run  through  the  truth  model  and  estimator.  Each  test  case 
produced  ten  data  sets  with  different  values  of  noise.  Table  12  lists  the  magnitudes  and 
individual  components  of  non-gravitational  acceleration  used  for  each  test  case.  Large 
accelerations  such  as  A  g  or  greater  usually  produce  noticeable  results;  therefore,  it  is  the 
smaller  accelerations  of  magnitude  cm/s"  which  are  of  particular  interest  in  this  research. 
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"able  12.  IN 

on-Gravitational  Acceleration  Present 

Case  1 

Case  2 

Case  3 

Case  4 

Case  5 

Case  6 

A0x  (km/s2) 

1.00E-04 

1.00E-03 

2.00E-04 

9.00E-05 

4.00E-06 

2.89E-06 

A0y  (km/s2) 

-2.00E-05 

-1.00E-04 

8.00E-05 

2.00E-04 

3.00E-06 

-5.56E-06 

A0z  (km/s2) 

-3.00E-05 

-3.00E-04 

3.00E-04 

7.00E-05 

1.00E-05 

2.285E-06 

Magnitude  (km/s2) 

1.063E-04 

1.0489E-03 

3.69E-04 

2.30E-04 

1.12E-05 

6.67E-06 

The  values  listed  in  Table  12  were  input  into  the  truth  model  as  the  A(l  component. 

Table  13  displays  the  estimated  state  for  Case  1  after  the  Monte  Carlo  analysis.  Of 
special  interest  is  the  highlighted  section.  These  values  represent  the  estimated  amount  of 
total  acceleration  present  for  the  initial  state  vector.  The  estimated  value  due  to  non- 
gravitational  acceleration  is  found  by  subtracting  the  known  gravitational  acceleration 
value  found  in  the  previous  chapter,  (Equation  (73)  ).  Table  14  lists  the  estimated  values 
for  the  non-gravitational  acceleration  components. 


Table  13.  Case  1  Monte  Carlo  Estimated  State 


Component 

Estimated  State 

x  (km) 

-6079.601784 

y  (km) 

1837.901274 

z  (km) 

-1596.595273 

X  (km/s) 

-2.959936964 

y  (km/s) 

-5.6501156 

Z  (km/s) 

4.820181948 

ax  (km/s2) 

0.008737164 

ax  (km/s2) 

-0.002630849 

a.  (km/s2) 

0.002243172 

ax  (km/s3) 

4.2860  IE-06 

aY  (km/s3) 

8.03442E-06 

d_(km/s3) 

-6.941 13E-06 
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Table  14.  Case  1  Estimated  Non-gravitational  Acceleration 


Case  1 

Ax  (km/s2) 

0.000100007 

Ay  (km/s2) 

-1.97844E-05 

Az  (km/s2) 

-3.20633E-05 

Listed  in  Table  15  are  the  values  for  the  variance  from  the  Monte  Carlo 
covariance  matrix.  The  standard  deviation  (a)  obtained  from  these  variances  is  also  listed 
in  Table  15,  as  well  as  the  estimated  minimum  and  maximum  values  for  non-gravitational 
acceleration,  which  were  determined  by  taking  the  estimated  value  from  Table  14  and 
both  subtracting  and  adding  the  standard  deviation,  respectively.  By  inspection,  it  is 
apparent  that  the  true  value  is  located  between  the  minimum  and  maximum  estimates.  So 
far,  the  truth  model  and  estimator  have  yielded  the  desired  results. 


Table  15.  Case  1  Monte  Carlo  Results 


Variable 

Variance 

Sigma  (o) 

Estimate  (min) 

Estimate  (max) 

True 

Ax  (km/s2) 

1.1 15E-10 

1.05594E-05 

8.94475E-05 

1.10566E-04 

1.00E-04 

Ay  (km/s2) 

1.0278E-1 1 

3.20593E-06 

-1.65785E-05 

-2.29904E-05 

-2.00E-05 

A0:  (km/s2) 

1.7407E-1 1 

4.17217E-06 

-2.78911E-05 

-3.62355E-05 

-3.00E-05 

The  process  used  to  obtain  the  results  in  Tables  13-15  for  Case  1  was  used  for  all 
subsequent  test  cases.  Table  16  displays  the  estimated  state  vectors  for  all  test  cases 
using  the  Monte  Carlo  analysis.  The  highlighted  rows  are  the  various  estimated 
components  of  a  with  gravitational  and  non-gravitational  acceleration  combined. 
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Table  16.  Monte  Carlo  Estimated  State  Vector 


Case  1 

Case  2 

Case  3 

Case  4 

Case  5 

Case  6 

r 

-6079.601784 

-6079.592763 

-6079.617408 

-6079.595845 

-6079.605393 

-6079.617422 

1837.901274 

1837.899984 

1837.922794 

1837.897105 

1837.910345 

1837.92279 

-1596.595273 

-1596.602284 

-1596.609607 

-1596.597832 

-1596.610173 

-1596.609597 

V 

-2.959936964 

-2.959752536 

-2.959699776 

-2.960101323 

-2.959403419 

-2.959698775 

-5.6501156 

-5.649968385 

-5.650207202 

-5.650044755 

-5.650208869 

-5.650207324 

4.820181948 

4.819868758 

4.819932365 

4.820190012 

4.819921725 

4.819932347 

a 

0.008737164 

0.009630106 

0.008843116 

0.008726374 

0.008641327 

0.008646013 

-0.002630849 

-0.002711875 

-0.00253619 

-0.002409998 

-0.002608829 

-0.002621751 

0.002243172 

0.001973669 

0.002576583 

0.002343208 

0.002287571 

0.002278865 

a 

4.28601E-06 

4.16182E-06 

4.19147E-06 

4.36836E-06 

4.03754E-06 

4.19097E-06 

8.03442E-06 

8.08126E-06 

8.09605E-06 

7.99999E-06 

8.07227E-06 

8.09606E-06 

-6.94113E-06 

-6.82963E-06 

-6.83764E-06 

-6.95477E-06 

-6.81189E-06 

-6.83767E-06 

Just  like  Case  1,  the  known  gravitational  acceleration  value  is  subtracted  from  a  to 
obtain  the  estimated  non-gravitational  acceleration  values  (Table  17).  Tables  18-21  list 
the  results  for  Cases  2-5. 


Table  17.  Estimated  Non-gravitational  Acceleration  for  all  Cases 


Case  1 

Case  2 

Case  3 

Case  4 

Case  5 

Case  6 

X 

0.000100007 

0.000992949 

0.00020596 

8.92169E-05 

4.17071E-06 

8.85662E-06 

X 

-1.97844E-05 

-0.00010081 

7.48753E-05 

0.000201067 

2.23559E-06 

-1.06865E-05 

X 

-3.20633E-05 

-0.000301566 

0.000301348 

6.79723E-05 

1.23354E-05 

3.62937E-06 

Tab 

e  18.  Case  2  IV 

onte  Carlo  Results 

Variable 

Variance 

Sigma  (a) 

Estimate 

Estimate 

True  Value 

X  (km/s2) 

6.3 165E- 1 1 

7.94764E-06 

9.85001E-04 

1.000897E-03 

1.00E-03 

A0v  (km/s2) 

3.71 16E-1 1 

6.09229E-06 

-1.06902E-04 

-9.47179E-05 

-1.00E-04 

X  (km/s2) 

2.1 172E-11 

4.6013E-06 

-3.06167E-04 

-2.96965  E-04 

-3.00E-04 
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Tab] 

le  19.  Case  3  IV 

ionte  Carlo  Results 

Variable 

Variance 

Sigma  (o) 

Estimate 

Estimate 

True  Value 

Ax  (km/s2) 

1.027E-10 

1.01341E-05 

1.95826E-04 

2.16094E-04 

2.00E-04 

Ay  (km/s2) 

9.9128E-11 

9.9563E-06 

6.4919E-05 

8.483 16E-05 

8.00E-05 

Az  (km/s2) 

4.188E-11 

6.47148E-06 

2.94876E-04 

3.07819E-04 

3.00E-04 

Tab] 

le  20.  Case  4  IV 

ionte  Carlo  Results 

Variable 

Variance 

Sigma  (o) 

Estimate 

Estimate 

True  Value 

Ax  (km/s2) 

9.908E-11 

9.95389E-06 

7.9263E-05 

9.91708E-05 

9.00E-05 

Ay  (km/s2) 

1.1389E-11 

3.37476E-06 

1.97692  E-04 

2.04442  E-04 

2.00E-04 

Az  (km/s2) 

1.7184E-11 

4.14536E-06 

6.38269E-05 

7.21 176E-05 

7.00E-05 

Tab] 

le  2 1 .  Case  5  IV 

ionte  Carlo  Results 

Variable 

Variance 

Sigma  (<t) 

Estimate 

Estimate 

True  Value 

Ax  (km/s2) 

4.3922E-1 1 

6.62737E-06 

-2.45665E-06 

1.07981E-05 

4.00E-06 

Ay  (km/s2) 

4.4168E-1 1 

6.6459E-06 

-4.41031E-06 

8.88149E-06 

3.00E-06 

Az  (km/s2) 

1.7455E-11 

4.17792E-06 

8.1575E-06 

1.65133E-05 

1.00E-05 

There  are  several  noteworthy  outcomes  obtained  from  Cases  1-5.  First  of  all,  the 
values  obtained  in  Table  17  for  Cases  1-5  are  extremely  close  to  the  true  non- 
gravitational  acceleration  values  input  into  the  truth  model.  It  is  also  important  to  note 
that  the  minimum  and  maximum  values  found  by  taking  into  account  the  standard 
deviation  are  also  quite  close  to  the  true  value.  The  results  from  these  cases  support  the 
validity  of  the  estimator. 

The  results  for  Case  6  are  quite  different  than  those  seen  in  the  previous  five 
cases.  The  non-gravitational  acceleration  values  in  Table  17  do  not  reflect  the  true  values 
whatsoever.  The  standard  deviation  is  also  the  highest  it  has  been  for  any  case.  In  fact, 
for  the  A0  component,  the  standard  deviation  has  a  greater  magnitude  than  the  estimated 

value.  When  the  standard  deviation  is  combined  with  the  values  in  Table  17,  there  is  no 
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clue  as  to  what  the  true  values  are  (Table  22).  The  estimator  is  incapable  of  estimating 
the  state  correctly  when  the  magnitude  of  non-gravitational  acceleration  is  that  small. 


Tab! 

le  22.  Case  6  IV 

ionte  Carlo  Results 

Variable 

Variance 

Sigma  (o) 

Estimate 

Estimate 

True  Value 

A0x  (km/s2) 

1.0287E-10 

1.01425E-05 

-1.28586E-06 

1.89991E-05 

2.89E-06 

A0y  (km/s2) 

9.942E-1 1 

9.97096E-06 

-2.06574E-05 

-7.15501E-07 

-5.56E-06 

A0:  (km/s2) 

4.2024E-1 1 

6.48259E-06 

-2.85322E-06 

1.01 12E-05 

2.285E-06 

Interestingly  enough,  the  results  from  Case  6  do  make  sense.  The  magnitude  of 
non-gravitational  acceleration  used  in  Case  6  is  the  smallest  magnitude  of  acceleration 
that  the  3rd  order  equations  of  motion  can  detect  given  the  conditions  outlined  in  the 
beginning  of  this  section.  This  is  true  if  the  constant  (C)  is  in  fact  zero.  The  constant, 
however,  is  more  than  likely  not  zero  but  some  other  value.  This  would  make  the 
magnitude  greater  than  0.667  cm/s  if  it  is  assumed  that  C  is  positive.  That  being  the 
case,  the  estimator  is  unable  to  accurately  estimate  the  state  at  such  small  magnitudes  due 
to  the  limitations  of  the  equations  of  motion  themselves. 


Real  World  Approach. 

The  above  analysis  would  make  it  seem  like  the  estimator  is  capable  of  detecting 
non-gravitational  acceleration  with  magnitude  as  small  as  1.12  cm/s";  unfortunately,  the 
method  used  to  obtain  the  above  results  does  not  reflect  reality.  In  the  real  world,  radar 
sites  do  not  have  N  number  of  data  sets  for  one  particular  satellite;  they  will  most  likely 
only  have  one  data  set  to  work  with.  This  being  the  case,  the  results  produced  by  the 
estimator  must  also  be  analyzed  using  a  real  world  approach.  Cases  1-5  will  be  analyzed 
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using  the  real  world  method.  Case  6  has  been  left  out  because  it  has  already  been  shown 
with  the  Monte  Carlo  analysis  that  the  results  are  invalid. 

Instead  of  averaging  the  ten  estimates  within  each  case  to  arrive  at  a  final 
estimate,  each  estimate  is  analyzed  as  an  individual  (N=l)  data  set.  For  each  of  these 
data  sets,  the  known  value  for  gravitational  acceleration  is  subtracted  from  a  to  obtain 
the  estimated  non-gravitational  acceleration  values.  Table  23  lists  these  estimated  non- 
gravitational  acceleration  values  for  five  of  the  data  sets  in  Case  1 .  Again,  in  reality  a 
radar  site  will  more  than  likely  only  have  one  data  set  to  work  with,  but  for  the  purposes 
of  this  research,  it  is  beneficial  to  see  if  the  correct  results  are  obtained  every  time  or  only 
once  in  a  while  by  chance.  After  comparing  the  estimated  values  to  their  true  values  in 
the  last  column  of  Table  23,  the  results  for  Case  1  appear  to  be  relatively  accurate; 
however,  various  components  in  some  of  the  data  sets  are  a  little  off.  For  example,  the 
estimated  value  of  the  A0y  component  for  set  1  is  closer  to  -  IE-05  km/s  and  in  set  2  and 

3  it  is  closer  to  -3E-05  km/s"  while  the  true  value  is  really  -2E-05  km/s". 

As  with  any  analysis,  it  is  important  to  take  into  account  the  standard  deviation  in 
order  to  get  a  true  understanding  of  the  accuracy  of  the  results.  The  standard  deviation 
values  used  in  this  section  were  obtained  from  the  MATLAB  filter  not  the  Monte  Carlo 
estimates.  Table  23  lists  the  ±  lo  and  ±  2a  estimates  for  the  various  non-gravitational 
acceleration  components  for  five  of  the  data  sets  in  Case  1 . 
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Table  23.  MATLAB  Case  1  Results  for  Individual  Data  Sets 


-2a 

-la 

Estimate 

la 

2a 

True 

Set  1 

Ar 

6.94671E-05 

7.8947E-05 

8.8428E-05 

9.79094E-05 

1.0739E-04 

1.00E-04 

Ay 

-2.7336E-05 

-2.0346E-05 

-1.3355E-05 

-6.36526E-06 

6.25089E-07 

-2.00E-05 

A, 

-4.0340E-05 

-3.4700E-05 

-2.9060E-05 

-2.34206E-05 

-1.7780E-05 

-3.00E-05 

Set  2 

A, 

8.27126E-05 

9.2193E-05 

1.0167E-04 

1.1 1 154E-04 

1.20635E-04 

1.00E-04 

Ay 

-4.21 18E-05 

-3.5127E-05 

-2.8136E-05 

-2.1 1452E-05 

-1.4154E-05 

-2.00E-05 

AS 

-4.2266E-05 

-3.6626E-05 

-3.0986E-05 

-2.5347E-05 

-1.9707E-05 

-3.00E-05 

Set  3 

\ 

9.09363E-05 

1.0041E-04 

1.0989E-04 

1.19379E-04 

1.28859E-04 

1.00E-04 

Ay 

-4. 1090E-05 

-3.4100E-05 

-2.7109E-05 

-2.01198E-05 

-1.3129E-05 

-2.00E-05 

A: 

-3.9096E-05 

-3.3456E-05 

-2.7816E-05 

-2.21768E-05 

-1.6536E-05 

-3.00E-05 

Set  4 

A, 

8.61503E-05 

9.5630E-05 

1.0511E-04 

1.14592E-04 

1.24072E-04 

1.00E-04 

Ay 

-3.4771E-05 

-2.7780E-05 

-2.0790E-05 

-1.37995E-05 

-6.8089E-06 

-2.00E-05 

A: 

-4.5425E-05 

-3.9785E-05 

-3.4146E-05 

-2.85064E-05 

-2.2866E-05 

-3.00E-05 

Set  5 

A, 

9.08761E-05 

1.0035E-04 

1.0983E-04 

1 . 1 93 1 8E-04 

1.28799E-04 

1.00E-04 

Ay 

-3.0967E-05 

-2.3977E-05 

-1.6986E-05 

-9.99657E-06 

-3.0062E-06 

-2.00E-05 

A: 

-3.6229E-05 

-3.0589E-05 

-2.495E-05 

-1.93102E-05 

-1.3670E-05 

-3.00E-05 

By  inspection,  there  are  several  components  where  the  estimate  does  not  contain 
the  true  value  within  ±  la.  These  components  have  been  highlighted  for  easier 
identification.  These  highlighted  components  happen  to  be  off  by  a  magnitude  of  only 
10'  km/s  or  10'  km/s  .  This  error  is  minimal  considering  the  true  values  are  of 
magnitude  10'4  km/s2  and  1 0‘"  km/s2. 

Of  the  15  components  listed  in  Table  23,  five  of  them  do  not  contain  their  true 
value  within  ±  la.  In  other  words,  roughly  33%  of  the  answers  are  not  within  ±  la.  In 
Gaussian  statistics  the  probability  that  the  answer  is  not  within  ±  la  is  32%,  therefore,  the 
distribution  of  the  results  in  Case  1  is  expected.  All  data  sets  in  Case  1  are  accurate 
within  ±  2a. 
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The  results  for  Cases  2  and  3  are  quite  similar  to  Case  1 .  The  estimated  values  in 
Table  24  for  Case  2  and  Table  25  for  Case  3  are  pretty  consistent  across  the  data  sets. 
Compared  to  their  true  values  in  the  final  column  of  Tables  24  and  25,  the  estimates 
appear  to  be  very  close.  In  fact,  Cases  2  and  3  appear  to  be  better  at  estimating  the  true 
accelerations  than  Case  1 .  This  may  be  because  Cases  2  and  3  have  larger  non- 
gravitational  acceleration  components,  which  makes  it  harder  for  their  values  to  get  lost 
in  the  noise. 


Table  24.  MATLAB  Case  2  Results  for  Individual  Data  Sets 


-2c 

-la 

Estimate 

la 

2a 

True 

Set  1 

Ax 

9.81200E-04 

9.90686E-04 

1.00017E-03 

1.00966E-03 

1.01914E-03 

1.00E-03 

Ay 

-1.0891E-04 

-1.0192E-04 

-9.4929E-05 

-8.7936E-05 

-8.0943E-05 

-1.00E-04 

A; 

-3.1595E-04 

-3.1031E-04 

-3.0467E-04 

-2.9903E-04 

-2.9340E-04 

-3.00E-04 

Set  2 

Ax 

9.85444E-04 

9.94930E-04 

1.00442E-03 

1.01390E-03 

1.02339E-03 

1.00E-03 

Ay 

-1.3704E-04 

-1.3005E-04 

-1.2305E-04 

-1. 1606E-04 

-1.0907E-04 

-1.00E-04 

A: 

-3.0151E-04 

-2.9587E-04 

-2.9024E-04 

-2.8460E-04 

-2.7896E-04 

-3.00E-04 

Set  3 

Ax 

9.87489E-04 

9.96975E-04 

1.00646E-03 

1.01595E-03 

1.02543E-03 

1.00E-03 

Ay 

-1.2624E-04 

-1.1925E-04 

-1.1225E-04 

-1.0526E-04 

-9.8271E-05 

-1.00E-04 

A; 

-3.01 15E-04 

-2.9551E-04 

-2.8987E-04 

-2.8424E-04 

-2.7860E-04 

-3.00E-04 

Set  4 

Ax 

9.9323  7E-04 

1.00272E-03 

1.01221E-03 

1.02169E-03 

1 .03 1 1 8E-03 

1.00E-03 

Ay 

-1.1915E-04 

-1.1216E-04 

-1.0517E-04 

-9.8176E-05 

-9.1183E-05 

-1.00E-04 

A: 

-3.0529E-04 

-2.9965E-04 

-2.9401E-04 

-2.8837E-04 

-2.8273E-04 

-3.00E-04 

Set  5 

Ax 

1.00255E-03 

1.01204E-03 

1.02152E-03 

1.03 101E-03 

1.04049E-03 

1.00E-03 

Ay 

-1.2343E-04 

-1. 1643E-04 

-1.0944E-04 

-1.0245E-04 

-9.5458E-05 

-1.00E-04 

A- 

-3.2194E-04 

-3.1630E-04 

-3.1066E-04 

-3.0502E-04 

-2.9938E-04 

-3.00E-04 

Tables  24  and  25  also  take  into  account  the  standard  deviation  for  Cases  2  and  3. 
As  seen  with  Case  1,  most  of  these  data  sets  have  a  component  where  the  estimated  value 
is  not  within  ±  la  of  the  true  value.  Again,  these  components  have  been  highlighted. 
These  components  are  off  by  a  magnitude  of  10'  km/s  or  10'  km/s  error.  Since  the 
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true  values  are  of  magnitude  I  O'3  km/s2  through  10’5  km/s2,  the  magnitude  of  error  does 
not  necessarily  destroy  the  entire  validity  of  the  results.  It  is  still  easy  to  ascertain  the 
magnitude  of  the  true  values. 


Table  25.  MATLAB  Case  3  Results  for  Individual  Data  Sets 


-2c 

-la 

Estimate 

la 

2a 

True 

Set  1 

A, 

1.84874E-04 

1.94345E-04 

2.03815E-04 

2.13286E-04 

2.22757E-04 

2.00E-04 

Ay 

7.09803E-05 

7.79598E-05 

8.49393E-05 

9.19188E-05 

9.88982E-05 

8.00E-05 

A; 

2.86076E-04 

2.9171  IE-04 

2.97346E-04 

3.02982E-04 

3.08617E-04 

3.00E-04 

Set  2 

A r 

1.62638E-04 

1.72109E-04 

1.81579E-04 

1.9105E-04 

2.0052E-04 

2.00E-04 

Ay 

6.83558E-05 

7.53357E-05 

8.23156E-05 

8.92955E-05 

9.62754E-05 

8.00E-05 

A: 

2.8615  IE-04 

2.91786E-04 

2.97421E-04 

3.03056E-04 

3.08691E-04 

3.00E-04 

Set  3 

A r 

1.75432E-04 

1.84902E-04 

1.94373E-04 

2.03843E-04 

2.13314E-04 

2.00E-04 

Ay 

6.40875E-05 

7.10672E-05 

7.80468E-05 

8.50265E-05 

9.20062E-05 

8.00E-05 

A: 

2.9456E-04 

3.00195E-04 

3.05829E-04 

3.1 1464E-04 

3.17099E-04 

3.00E-04 

Set  4 

A, 

1.97613E-04 

2.07083E-04 

2.16553E-04 

2.26023E-04 

2.35493E-04 

2.00E-04 

Ay 

6.4355  IE-05 

7.13347E-05 

7.83 143E-05 

8.52939E-05 

9.22735E-05 

8.00E-05 

A: 

2.88561E-04 

2.9983  IE-04 

2.94196E-04 

2.82926E-04 

3.05466E-04 

3.00E-04 

Set  5 

A, 

1.7555E-04 

1.94491E-04 

1.8502E-04 

1.66079E-04 

2.03961E-04 

2.00E-04 

Ay 

7.36103E-05 

8.75702E-05 

8.059E-05 

6.66303E-05 

9.45502E-05 

8.00E-05 

A: 

2.91338E-04 

3.02608E-04 

2.9697E-04 

2.85703E-04 

3.08243E-04 

3.00E-04 

As  seen  with  Cases  1-3,  the  estimates  in  Table  26  for  Case  4  are  consistent  across 
the  data  sets  and  are  close  to  portraying  the  true  values.  In  fact,  the  results  in  Table  26 
seem  to  be  the  best.  Data  set  5  is  the  only  set  that  has  a  component  not  within  ±  la  of  the 
true  value. 
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Table  26.  MATLAB  Case  4  Results  for  Individual  Data  Sets 


-2a 

-lo 

Estimate 

lo 

la 

True 

Set  1 

Ax 

8.01006E-05 

8.95745E-05 

9.90484E-05 

1.08522E-04 

1.17996E-04 

9.00E-05 

Ay 

1.81295E-04 

1.88282E-04 

1.95268E-04 

2.02255E-04 

2.0924  IE-04 

2.00E-04 

A: 

5.89685E-05 

6.46025E-05 

7.02365E-05 

7.58705E-05 

8.15045E-05 

7.00E-05 

Set  2 

Ax 

7.24923E-05 

8.19667E-05 

9.1441E-05 

1.00915E-04 

1.1039E-04 

9.00E-05 

Ay 

1.89723E-04 

1.96709E-04 

2.03695E-04 

2.10681E-04 

2.1766E-04 

2.00E-04 

A; 

5.33809E-05 

5.90152E-05 

6.46495E-05 

7.02837E-05 

7.5918E-05 

7.00E-05 

Set  3 

A, 

7.45168E-05 

8.39912E-05 

9.34655E-05 

1.0294E-04 

1.12414E-04 

9.00E-05 

Ay 

1.85354E-04 

1.9234E-04 

1.99326E-04 

2.063 13E-04 

2.13299E-04 

2.00E-04 

A , 

6.01894E-05 

6.58237E-05 

7.1458E-05 

7.70922E-05 

8.27265E-05 

7.00E-05 

Set  4 

A, 

6.78759E-05 

7.73498E-05 

8.68237E-05 

9.62976E-05 

1.05772E-04 

9.00E-05 

Ay 

1.8821E-04 

1.95196E-04 

2.02183E-04 

2.09169E-04 

2.16156E-04 

2.00E-04 

A: 

5.72881E-05 

6.29223E-05 

6.85566E-05 

7.41909E-05 

7.98251E-05 

7.00E-05 

Set  5 

Ax 

7.11074E-05 

8.05815E-05 

9.00557E-05 

9.95298E-05 

1.09004E-04 

9.00E-05 

Ay 

1.9201E-04 

1.98996E-04 

2.05982E-04 

2.12968E-04 

2.19954E-04 

2.00E-04 

A: 

5.17785E-05 

5.74129E-05 

6.30472E-05 

6.86816E-05 

7.43 159E-05 

7.00E-05 

In  the  previous  section,  Case  5  contained  the  smallest  magnitude  of  non- 
gravitational  acceleration  that  the  estimator  could  detect  using  the  Monte  Carlo  method. 
The  results  in  Table  27,  however,  show  that  the  estimator  cannot  find  the  estimate 
reliably.  Not  a  single  set  has  estimated  values  that  are  close  to  the  true  non-gravitational 
acceleration  values.  The  true  value  may  have  been  obtainable  after  taking  an  average 
using  Monte  Carlo  analysis;  but  as  separate  individual  estimates  the  true  values  are 
unclear. 
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ATLAB  Case  5  Results  for  Individual  Data  Sets 


Table  27.  M 


-2a 

-lc 

Estimate 

la 

2  a 

True 

Set  1 

Ar 

-3.2597E-06 

6.2199E-06 

1.5699E-05 

2.5179E-05 

3.4658E-05 

4.00E-06 

o 

-1.4599E-05 

-7.6089E-06 

-6.189E-07 

6.371  IE-06 

1.3361E-05 

3.00E-06 

A: 

-4.2422E-06 

1.3977E-06 

7.0376E-06 

1.2677E-05 

1.8317E-05 

1.00E-05 

Set  2 

Ax 

-2.4232E-05 

-1.4752E-05 

-5.2729E-06 

4.2068E-06 

1.3686E-05 

4.00E-06 

o 

9.6891E-07 

7.9588E-06 

1.4948E-05 

2.1938E-05 

2.8928E-05 

3.00E-06 

A: 

-4.1779E-06 

1.4619E-06 

7.1017E-06 

1.2741E-05 

1 .838  IE-05 

1.00E-05 

Set  3 

Ar 

-2.381E-05 

-1.4330E-05 

-4.8512E-06 

4.6280E-06 

1.4107E-05 

4.00E-06 

Ay 

-9.0425E-06 

-2.0522E-06 

4.9380E-06 

1.1928E-05 

1.8918E-05 

3.00E-06 

A: 

-6.0754E-06 

-4.3569E-07 

5.2040E-06 

1.0843E-05 

1.6483E-05 

1.00E-05 

Set  4 

A, 

-2.0674E-05 

-1.1 195E-05 

-1.7156E-06 

7.7639E-06 

1.7243E-05 

4.00E-06 

o 

-2.4414E-05 

-1.7424E-05 

-1.0434E-05 

-3.4442E-06 

3.5457E-06 

3.00E-06 

A: 

-1.7366E-06 

3.9030E-06 

9.5427E-06 

1.5182E-05 

2.0822E-05 

1.00E-05 

Set  5 

Ax 

-2.0173E-05 

-1.0694E-05 

-1.2148E-06 

8.2646E-06 

1.7744E-05 

4.00E-06 

Ay 

-1.3768E-05 

-6.7779E-06 

2.1217E-07 

7.2023E-06 

1.4192E-05 

3.00E-06 

A , 

-3.4714E-06 

2.1676E-06 

7.8068E-06 

1.3446E-05 

1.9085E-05 

1.00E-05 

Table  27  also  lists  the  values  for  ±  la  and  ±  2a.  At  first  glance,  it  may  seem  that 
since  there  are  only  three  highlighted  components  the  results  must  be  fairly  good,  but  this 
is  not  the  case.  The  amount  of  error  between  these  highlighted  components  to  their  true 
values  is  of  magnitude  10'  km/s  .  This  also  happens  to  be  the  same  magnitude  of  the 
true  values  themselves.  Ultimately,  what  this  means  is  that  the  estimator  has  no  idea 
what  the  true  values  really  are  and  does  a  poor  job  estimating  at  such  a  small  magnitude. 

Depending  on  the  actual  use  of  the  estimator,  results  such  as  these  may  still  be 
useful.  The  estimator  is  still  indicating  that  there  is  non-gravitational  acceleration  present 
even  if  the  estimate  is  not  statistically  accurate.  Sometimes  it  is  more  important  to 
assume  an  object  has  non-gravitational  acceleration  present  when  it  does  not  than  to 
assume  an  object  does  not  have  non-gravitational  acceleration  when  it  does.  Further 
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analysis  can  always  be  conducted  on  the  objects  that  are  assumed  to  contain  non- 
gravitational  acceleration. 

After  careful  analysis  of  the  previous  five  cases,  the  following  conclusions  can  be 
drawn.  The  uncertainty  in  the  results  for  Case  5  is  too  high  for  the  estimator  to  be  of 
much  use  at  such  a  small  magnitude  of  non-gravitational  acceleration.  The  results  seen  in 
Cases  1-4  are  much  more  accurate.  Of  these  four  cases,  Case  1  had  the  smallest 
magnitude  at  10.63  cm/s  .  Based  on  the  above  analysis  and  results,  it  is  fairly  safe  to  say 
that  the  estimator  can  detect  non-gravitational  acceleration  down  to  a  magnitude  of  10.63 
cm/s“.  Non-gravitational  acceleration  that  only  has  components  of  magnitude  10'  km/s 
and  10'“  km/s"  seem  to  get  lost  in  both  the  noise  and  the  estimation  capability  of  the 
estimator.  In  these  cases,  the  estimator  is  able  to  detect  the  non-gravitational  acceleration 
but  is  unable  to  give  a  truly  decent  estimate  as  to  its  true  value. 

The  inability  to  measure  these  small  magnitudes  may  be  caused  by  a  combination 
of  factors.  First  and  foremost,  the  general  equations  used  to  model  the  motion  of  the 
space  object  already  limits  the  detection  level  of  the  results.  Also,  it  is  important  to 
realize  that  there  are  going  to  be  inaccuracies  caused  by  MATLAB  itself.  MATLAB,  as 
with  many  other  computer  programs,  will  ultimately  truncate  numbers  as  they  are  being 
run  through  the  computer  code.  This  truncation  may  happen  in  the  10th  decimal  place  or 
even  higher  and  may  seem  insignificant,  but  can  produce  very  noticeable  results 
especially  when  working  with  such  small  acceleration  values. 
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Investigative  Questions  Answered 

Several  questions  were  posed  in  Chapter  I  concerning  the  development  of  the 
truth  model  and  the  detection  capability  of  the  estimator.  The  first  question,  concerning 
how  much  non-gravitational  acceleration  the  estimator  would  be  able  to  detect  is 
answered  after  analyzing  numerous  test  cases.  It  was  shown  that  the  estimator  should  be 
able  to  detect  as  little  as  10.63  cm/s  in  magnitude.  The  second  question  dealt  with 
Galileo’s  projectile  equations  and  whether  or  not  they  would  be  accurate  enough  to 
model  an  object  in  space.  It  was  shown  that  these  equations  were  in  fact  too  general  to 
obtain  sufficient  results;  however,  a  Taylor  series  expansion  to  the  3rd  order  achieved 
success.  The  third  question  was  concerning  the  geopotential.  Multiple  terms  such  as  J2, 
J4  or  two-body  effects  can  go  into  determining  the  geopotential.  For  this  research,  only 
two-body  and  Jo  effects  comprised  the  terms  of  the  geopotential.  These  two  terms  proved 
quite  effective  at  modeling  the  geopotential  and  its  gradient.  Proof  can  be  seen  in 
Equations  (73)  and  (74)  where  the  MATLAB  and  STK  acceleration  values  at  epoch  are 
nearly  identical. 

Summary 

Various  methods  were  used  to  analyze  the  simulated  data  to  obtain  valid  results. 
An  STK  simulation  was  used  as  a  baseline  model  to  verify  the  accuracy  of  the  equations 
of  motion.  The  2nd,  3rd,  and  4th  order  Taylor  series  approximations  were  all  tested  to  see 
which  would  yield  better  results  at  modeling  an  orbiting  space  object.  The  3ld  order 
equations  of  motion  were  shown  to  have  the  required  accuracy  necessary  for  detecting 
non-gravitational  acceleration.  Numerous  test  cases  (both  with  and  without  non- 
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gravitational  acceleration  present)  were  tested  with  the  estimator.  The  Monte  Carlo 
method  was  used  to  verify  the  statistics  of  the  results.  A  real  world  approach  was  used  to 
analyze  the  various  scenarios  to  determine  the  smallest  magnitude  that  could  be 
determined  with  certainty.  Using  the  least  squares  estimator,  accelerations  with  a 
magnitude  as  small  as  10.63  cm/s2  were  detectable  with  statistical  accuracy.  At 
magnitudes  of  acceleration  smaller  than  10.63  cm/s  ,  confidence  and  validity  of  the 
results  declines.  Nearly  all  estimates  were  within  2a;  however,  this  is  not  beneficial  in 
the  test  cases  where  the  magnitude  of  standard  deviation  is  the  same  as  or  greater  than  the 
magnitude  of  the  estimated  values. 
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V.  Conclusions  and  Recommendations 


Chapter  Overview 

The  purpose  of  this  chapter  is  to  conclude  the  findings  and  discuss  the 
significance  of  this  research.  Suggestions  as  to  what  actions  should  be  taken  as  a  result 
of  the  findings  will  be  explored.  Various  aspects  of  further  research  will  also  be 
addressed. 

Conclusions  of  Research 

It  has  been  shown  that  a  3rd  order  Taylor  series  expansion  can  adequately  model 
objects  in  LEO  for  a  short  period  of  time.  These  equations  of  motion,  with  the  use  of  a 
truth  model  and  linear  least  squares  estimator,  can  detect  constant  non-gravitational 
acceleration  down  to  roughly  10.63  cm/s"  in  magnitude  with  statistical  accuracy.  If  the 
amount  of  non-gravitational  acceleration  is  smaller  than  the  10.63  cm/s"  limit,  the 
estimator  produces  poor  results.  At  the  very  least,  the  estimator  is  capable  of  detecting 
non-gravitational  acceleration  but  may  not  necessarily  be  able  to  determine  the  exact 
magnitude  of  its  components  within  la. 

Significance  of  Research 

There  are  several  significant  aspects  of  this  research.  First  of  all,  the  above 
research  has  shown  that  the  motion  of  space  objects  can  be  modeled  using  linear 
dynamics  for  a  short  time  span.  Adding  non-gravitational  acceleration  to  the  dynamics, 
(as  seen  with  objects  such  as  tethered  systems,  maneuvering  satellites,  and  thrusting 
ballistic  missiles)  the  effects  of  non-Keplerian  motion  can  now  be  tracked  using  linear 
equations.  Even  the  data  in  range,  azimuth,  and  elevation  fonnat  can  be  exchanged  for 
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linear  data  in  the  fonn  of  position  and  velocity.  The  use  of  a  linear  model  enables 
estimates  of  the  initial  state  to  be  calculated  at  the  click  of  a  button.  This  is  much  more 
favorable  than  that  of  iterative  methods  seen  in  nonlinear  systems. 

There  are  currently  a  plethora  of  elaborate  methods  available  for  identifying  and 
tracking  objects  that  are  solely  tethered  or  solely  ballistic  missiles;  however,  there  are  not 
many  methods  that  follow  a  more  general  approach  and  distinguish  simply  between 
Keplerian  and  non-Keplerian  objects.  This  general  approach  allows  for  the  quicker 
detection  of  non-Keplerian  objects.  Then,  if  necessary,  a  more  in-depth  analysis  can  be 
used  to  identify  whether  the  object  is  a  tethered  system,  an  object  on  re-entry,  or  a 
maneuvering  satellite. 

Recommendations  for  Action 

It  is  recommended  that  radar  sites  utilize  this  computer  model  in  order  to  filter  out 
space  objects  that  contain  a  certain  level  of  non-gravitational  acceleration.  This  computer 
model  is  not  able  to  track  objects  for  an  extended  period  of  time,  nor  is  it  able  to 
distinguish  between  a  tethered  system  and  an  object  on  re-entry;  however,  it  is  quite 
useful  as  a  filter.  The  model  is  rather  effective  at  determining  the  state  of  an  object  with 
or  without  non-gravitational  acceleration  in  just  seconds.  This  is  very  important  when 
time  is  of  the  essence.  Depending  on  the  magnitude  of  the  non-gravitational  acceleration 
present,  other  identification  and  tracking  methods  can  be  used  to  provide  further  insight 
into  the  object  of  interest. 
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Recommendations  for  Future  Research 


There  are  various  avenues  of  future  research  dealing  with  the  current  research 
and,  in  general,  the  area  of  tracking  objects  with  non-Keplerian  motion.  The  above 
research  was  tested  and  validated  using  models  and  simulated  data.  It  would  be  quite 
beneficial  to  test  the  estimator  using  real-world  data.  Perhaps  testing  the  estimator  with 
data  from  satellites  that  have  on-board  accelerometers  could  be  advantageous.  This  way, 
there  is  still  some  indication  as  to  the  true  value  of  non-gravitational  acceleration  present. 
Another  recommendation  would  be  to  develop  a  model  that  can  detect  variable  non- 
gravitational  acceleration.  The  current  model  only  considers  constant  non-gravitational 
acceleration. 

Another  idea  to  consider  is  to  detennine  the  initial  state  vector  in  the  space 
object’s  coordinate  frame  instead  of  working  in  the  inertial  IJK  coordinate  frame.  This 
approach  would  provide  some  very  important  insights.  For  instance,  a  lot  can  be  learned 
about  a  tethered  system  simply  given  the  various  components  of  non-gravitational 
acceleration  in  the  satellite’s  coordinate  frame.  Non-gravitational  acceleration  in  a 
tethered  system  can  only  be  located  in  the  tangential  or  radial  directions.  If  the 
magnitude  in  both  of  these  directions  is  zero,  then  the  object  is  not  tethered.  If  the 
tangential  is  zero  but  the  radial  is  not,  then  it  can  be  deduced  that  the  system  is  oriented 
vertically.  If  the  radial  is  zero  and  the  tangential  is  not,  then  the  system  is  horizontally 
oriented.  If  neither  component  is  zero,  then  the  tether  orientation  angle,  magnitude  of 
tether  force,  apparent  gravitational  parameter,  and  distance  to  the  center  of  mass  can  be 
determined  (Cicci  and  others,  2001a:  316-317).  Unfortunately,  these  statements  are  true 
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only  if  the  tethered  system  is  moving  in  the  plane.  Out  of  plane  motion  requires  a 
different  approach. 

Summary 

Objects  that  contain  non-gravitational  acceleration;  such  as  tethered  systems, 
thrusting  ballistic  missiles  and  maneuvering  satellites;  can  be  modeled  using  a  3rd  order 
Taylor  series  expansion.  The  combination  of  these  dynamics  with  a  linear  least  squares 
estimator  provides  the  ability  to  accurately  estimate  the  initial  state  and  covariance  of  an 
object  in  space.  The  sum  of  gravitational  and  non-gravitational  acceleration  acting  on  the 
object  is  part  of  the  state  vector.  A  non-zero  magnitude  for  non-gravitational  acceleration 
means  that  the  space  object  is  following  a  non-Keplerian  orbit.  The  current  model  used 
for  this  research  is  capable  of  detecting  non-gravitational  acceleration  as  small  as  10.63 
cm/s"  in  magnitude. 

The  identification  and  tracking  of  space  objects  will  remain  an  important  pursuit 
for  years  to  come.  As  technology  increases,  the  dynamics  of  a  space  object  will  more 
than  likely  not  follow  a  regular  Keplerian  orbit  due  to  the  desire  for  increased 
maneuverability  for  the  protection  against  potential  space  weapons.  Objects  such  as 
tethered  satellite  systems,  maneuvering  satellites  and  thrusting  ballistic  missiles  already 
follow  non-Keplerian  orbits.  Filters,  such  as  the  one  produced  in  this  research,  will 
become  a  necessity. 
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Appendix  A:  Equations  of  Motion 


Second  Order 

r  =  7i>  +  7oA^  +  ^  (A  +  go  ^ 
v  =  v0  +  (A0  +  g0)At 

a=(A+go) 

Third  Order 

1  -  -\X. 

r  =  r0  +  v0At  +  -(A0+g0)Ar  +-(A0  +  g0)At 

2  6 

v  =  v0  +  (4,  +  g0)At  +  +  |0)A^2 

«  =  (4)  +  fo)  +  (4>+ifo)Af 

a  =  (A0  +  g0 ) 

Fourth  Order 


r  =  r0  +  v0A?  +  ]-  ( A0  +  g0)At 2  +  ^(A0  +  g0)AF  +^~(A  +  go)At* 
2  6  24 

V  =  ^0  +  (A)  +  go  )A?  +  —  (A)  +  go  )' Ar  +  -  (A  +  go  )Af3 
2  6 

«  =  (A  +  g0 ) +  (A  +  go  )At + — (A  +  go  )■ Ar 
^  =  (4)+go)  +  (4)+go)A? 

«  =  (^o+#o) 
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Appendix  B:  MATLAB  Truth  Model 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%% 

%%  Author:  2Lt  Sandra  Rashash,  USAF,  17  May  07 

%% 

%% 

%%  Outputs: 


%% 

el  matrix 

-matrix  of  elevation  values 

rad 

%% 

az  matrix 

-matrix  of  azimuth  values 

rad 

%% 

rho  matrix 

-matrix  of  range  values 

km 

%% 

sitlat 

-radar  site  latitude 

rad 

%% 

sitlon 

-radar  site  longitude 

rad 

%% 

sitalt 

-radar  site  altitude 

km 

%% 

j  ulian  matrix 

-matrix  of  julian  days 

%% 

range  noise 

-range  noise 

km 

%% 

az  noise 

-azimuth  noise 

rad 

%% 

elnoise 

-elevation  noise 

rad 

%% 

%%  Locals: 

%% 

r  init 

-initial  position  at  epoch  in  IJK 

km 

%% 

X 

-T  component  of  position 

km 

%% 

y 

-T  component  of  position 

km 

%% 

z 

-'K'  component  of  position 

km 

%% 

v  init 

-initial  velocity  at  epoch 

km/s 

%% 

A 

-extra  acceleration  input 

km/sA2 

%% 

t 

-counter  for  time 

%% 

jd 

-time  in  julian  days 

%% 

gst 

-greenwhich  sidereal  time 

rad 

%% 

local  sideral  time 

-local  sidereal  time 

rad 

%% 

1st 

-matrix  of  local  sidereal  time 

rad 

%% 

time 

-time  since  epoch 

sec 

%% 

timematrix 

-matrix  of  times  since  epoch 

sec 

%% 

geopotential 

-geopotential 

mA2/sA2 

%% 

ax 

-'x'  component  of  acceleration 

km/sA2 

%% 

ay 

-'y'  component  of  acceleration 

km/sA2 

%% 

az 

-'z'  component  of  acceleration 

km/sA2 

%% 

g 

-local  gravity 

km/sA2 

%% 

r 

-new  position  vector 

km 

%% 

V 

-new  velocity  vector 

km/s 

%% 

i 

-index  variable 

%% 

vmatrix 

-matrix  of  velocities 

km/s 

%% 

x  site 

-station  coordinate  for  Ellipsoidal  Earth 

km 

%% 

z  site 

-station  coordinate  for  Ellipsoidal  Earth 

km 

%% 

Rsite 

-position  vector  of  radar  site 

km 

%% 

rot  ijk  sez 

-rotation  matrix  from  IJK  to  SEZ 

rad 

%% 

range  sez 

-range  vector  in  SEZ 

km 

%% 

rho  s 

-'s'  component  of  range 

km 

%% 

rho  e 

-'e'  component  of  range 

km 

%% 

rho  z 

-'z'  component  of  range 

km 

%% 

range  sez  matrix 

-matrix  of  range  values 

km 
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%% 

rho 

-range 

km 

%% 

az 

-azimuth 

rad 

%% 

el 

-elevation 

rad 

%% 

noise 

-matrix  of  noise  values 

%% 

ww 

-index  variable 

%% 

%% 

%%  Note:  An  extra  blank  line  is  present  at  the  end  of  the  output  file  that  must  be  deleted  before  sending 
the  data  to  the  estimator. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


WGS84Data 
format  long  g 

global  MU  J2  RE  TwoPI  EEarth  Rad 
radardata  1  =fopen('radardata  1  .dat',  'wt') ; 

%  site  data 

sitlat=  Rad*(-7.91); 
sitlon=Rad*(-14.4); 
sitalt=56.1*10A-3; 
range_noise=  101.7*10A-3; 
el_noise=Rad*(0.0283); 
az_noise=Rad*  (0 .0248); 

%  initial  conditions 
rinit  =[-6079. 6;  1837. 9;- 1596. 6]; 

x  =  r_init(l,l); 
y  =  r_init(2,l); 
z  =  r_init(3,l); 

v_init=[-2.96;-5.65;4.82]; 

u  =  v _init(l,l); 
v  =  v_init(2,l); 

w  =  v_init(3,l); 


A=[0;0;0]; 

A_dot=[0;0;0]; 


%  creates  a  5  minute  matrix  of  j  ulian  days  for  a  given  start  time 
julian_matrix(l,:)  =  JulianDay(2007,9, 13, 12,0,0); 

for  t=  1:299 

jd  =  JulianDay(2007,9,13,12,0,t); 

j  ulian_matrix(t+ 1 , :  )=jd; 
end 
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%  detemiines  matrix  for  the  local  sidereal  time 
for  index=l:300 


%  calculates  greenwich  sidereal  time 
gst  =  GSTime  (julian_matrix(index,:)); 

%  calculates  local  sidereal  time 
localsiderealtime  =  gst  +  sitlon; 

%  calls  revcheck  to  obtain  an  1st  less  than  2pi 

local  sidereal  time  =  rcvchecktlocal  sidcrcal  tirnc,  TwoPI); 

lst(index,:)=local_sidereal_time; 


%  converts  julian  days  into  time  since  epoch  observation 
time  =  (julian_matrix(index,:)-julian_matrix(151,l))*86400; 

%  creates  matrix  of  observation  times 
timematrix(  index,:)  =time; 

end 


%%  Determination  of  the  gravity  vector 
%  geopotential  taking  into  account  J2  and  2-body 
%  geopotential=(- 

MU/sqrt(xA2+yA2+zA2))+((MU*REA2*J2)/(2*(xA2+yA2+zA2)A(3/2)))*(3*zA2/(xA2+yA2+zA2)-l); 

% 

%  %  the  negative  gradient  of  the  geopotential 

% 

%  ax=-MU/(xA2+yA2+zA2)A(3/2)*x+15/2*zA2*MU*REA2*J2/(xA2+yA2+zA2)A(7/2)*x- 
3/2*MU*REA2*J2/(xA2+yA2+zA2)A(5/2)*x; 

% 

%  ay=-MU/(xA2+yA2+zA2)A(3/2)*y+15/2*zA2*MU*REA2*J2/(xA2+yA2+zA2)A(7/2)*y- 
3/2*MU*REA2*J2/(xA2+yA2+zA2)A(5/2)*y; 

% 

%  az=-MU/(xA2+yA2+zA2)A(3/2)*z- 

9/2*z*MU*REA2*J2/(xA2+yA2+zA2)A(5/2)+15/2*zA3*MU*REA2*J2/(xA2+yA2+zA2)A(7/2); 

% 

%  g=[ax;ay;az] 


%  2-body  acceleration  and  its  first  and  second  derivatives 
a_2body=  -MU*r_init*(mag(r_init)A-3); 

a_2body_dot=  3*MU*r_init*dot(r_init,v_init)*(mag(r_init)A-5)  -  MU*v_init*(mag(r_init)A-3); 
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%  J2  acceleration  and  its  first  and  second  derivatives 

aJ2=  [15/2*zA2*MU*REA2*J2/(xA2+yA2+zA2)A(7/2)*x-3/2*MU*REA2*J2/(xA2+yA2+zA2)A(5/2)*x;... 
15/2*zA2*MU*REA2*J2/(xA2+yA2+zA2)A(7/2)*y-3/2*MU*REA2*J2/(xA2+yA2+zA2)A(5/2)*y;... 
-9/2*z*MU*REA2*J2/(xA2+yA2+zA2)A(5/2)+15/2*zA3*MU*REA2*J2/(xA2+yA2+zA2)A(7/2)]; 

aJ2_dot=  -.5*MU*REA2*J2*[(-15*x*dot(r_init,v_init)*(mag(r_init)A-7)+3*u*(mag(r_init)A-5)... 
+105*x*zA2*dot(r_init,v_init)*(mag(r_init)A-9)-30*x*z*w*(mag(r_init)A-7)-15*u*zA2*(mag(r_init)A- 

7));...  ~  ' . “ 

(-15*y*dot(r_init,v_init)*(mag(r_init)A-7)+3*v*(mag(r_init)A-5)+105*y*zA2*dot(r_init,v_init)*... 

(mag(r_init)A-9)-30*y*z*w*(mag(r_init)A-7)-15*v*zA2*(mag(r_init)A-7));... 

(-45*z*dot(r_init,v_init)*(mag(r_init)A-7)+9*w*(mag(r_init)A- 

5)+105*zA3*dot(r_init,v_init)*(mag(r_init)A-9)-... 

45*zA2*w*(mag(r_init)A-7))]; 


%  total  acceleration  and  derivatives 
g=a_2body+aJ2 
g_dot=a_2body_dot+  aJ2_dot 


%%%  Equations  of  Motion  for  a  5  minute  radar  track  using  Taylor  series 
for  w=l:300 

t  =timematrix(w,:); 

%  finds  position 

r  =  r_init+  v_init*t  +((A+g)/2)*tA2+((A_dot+g_dot)/6)*tA3; 

R_ij  k_matrix(  w, :  )=r; 

%  finds  velocity 

v  =  v_init+(A+g)*t+((A_dot+g_dot)/2)*tA2; 
vmatrix(w,:)=v; 

%  finds  acceleration 

a  =  (A+g)  +  (A_dot+g_dot)*t; 

amatrix(w,:)=a; 

%  finds  1st  derivative  of  accel 
a_dot=(A_dot+g_dot) ; 
a_dot_matrix(w,:)=a_dot; 

end 


%  calculates  x  in  the  R_site  equation 

x_site  =  (RE/sqrt(l-EEarthA2*sin(sitlat)A2)+sitalt)*cos(sitlat); 

%  calculates  z  in  the  R  site  equation 

z_site  =  ((RE*(l-EEarthA2))/sqrt(l-EEarthA2*sin(sitlat)A2)+sitalt)*sin(sitlat); 


%  Determines  range,  azimuth,  and  elevation  matrices 
for  i=l:300 
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%  calculates  Rsite 

Rsite  =  [x_site*cos(lst(i));  x_site*sin(lst(i));  zsite]; 

%  rotation  matrix  from  the  IJK  to  SEZ  coord,  frame 

rot_ijk_sez  =  [sin(sitlat)*cos(lst(i))  sin(sitlat)*sin(lst(i))  -cos(sitlat);... 

-sin(lst(i))  cos(lst(i))  0;  cos(sitlat)*cos(lst(i))  cos(sitlat)*sin(lst(i))... 
sin(sitlat)]; 

range_sez=rot_ijk_sez*(R_ijk_matrix(i,:)'-R_site); 

range_sez_matrix(i,:)=range_sez; 

rhos  =range_sez_matrix(i,  1); 
rho_e  =range_sez_matrix(i,2); 
rhoz  =range_sez_matrix(i,3); 

%  determines  range 
rho=sqrt(rho_sA2+rho_eA2+rho_zA2); 
rho_matrix(i, :  )=rho ; 

%  determines  azimuth 
az  =  atan2(rho_e,-rho_s); 
if  az<0 

az=az+TwoPI; 

end 

az_matrix(i,  :)=az; 

%  determines  elevation 
el  =  asin(rho_z/rho); 
el_matrix(i,:)=el; 
end 

%  adds  noise  to  the  range, azimuth,  elevation  data  using  gaussian  random  number  generator 

noisel=randn(300,l)*.1017; 

noise2=randn(300,l)*. 000433; 

noise3=randn(300,l)*. 000494; 

%outputs  data  to  a  file 
for  ww=l:300 

fprintf(radardatal,'%22.15g',sitlat,sitlon,sitalt,julian_matrix(ww),rho_matrix(ww)+noisel(ww),az_matrix( 
ww)+noise2(ww),el_matrix(ww)+noise3(ww),range_noise,az_noise,el_noise); 
fprintf(radardata  1  ,'\n'); 

end 

%  %%  data  with  no  noise 
%  for  ww=l:300 

% 

fprintf(radardatal,'%22.15g',sitlat,sitlon,sitalt,julian_matrix(ww),rho_matrix(ww),az_matrix(ww),el_matrix 

(ww),range_noise,az_noise,el_noise); 

%  fprintf(radardatal,'\n'); 

% 

%  end 

fclose(radardata  1 ); 
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Appendix  C:  MATLAB  Estimator 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%% 

%% 


%%  Author:  2Lt  Sandra  Rashash,  USAF,  17  May  07 

%% 

%% 

%%  Locals: 

km 
sec 
rad 
rad 
rad 
km 
sec 

rad 
km 


%% 

phi  matrix 

-state  transition  matrix 

%% 

phi  multi  array 

-multi-dimensional  array  of  all  phi  matrices 

%% 

w 

-index  variable 

%% 

s 

-index  variable 

%% 

e 

-index  variable 

%% 

V 

-index  variable 

%% 

f 

-index  variable 

%% 

H 

-matrix  of  the  linearized  observation  relation 

%% 

T 

-observation  matrix 

%% 

%% 

%% 

T  matrix 

K 

K  array 

-multi-dimensional  array  of  all  T 

%% 

J 

-jacobian  matrix 

%% 

J  array 

-multi-dimensional  array  of  jacobians 

%% 

Qold 

-instrumental  covariance  for  rho,az,el 

%% 

Qnew 

-instrumental  covariance  for  position 

%% 

Q  array 

-multi-dimensional  array  for  all  Q  new 

%% 

c  variance 

-matrix  used  for  running  sums 

%% 

covariance 

-covariance  matrix 

%% 

s  vector 

-matrix  used  for  running  sums 

%% 

S 

-individual  state  vector  input  into  summation 

%% 

P 

-individual  covariance  input  into  summation 

%% 

%% 

state  vector 

-estimated  state  vector 

%% 

R  ijk  matrix 

-matrix  of  position  in  IJK  coords 

%% 

timematrix 

-matrix  of  observation  times 

%% 

1st 

-matrix  of  local  sidereal  time 

%% 

el  matrix 

-matrix  of  elevation 

%% 

az  matrix 

-matrix  of  azimuth 

%% 

rho  matrix 

-matrix  of  range 

%% 

t 

-time  variable 

%% 

jdtime 

-matrix  of  julian  days 

%% 

local  sideral  time 

-local  sidereal  time 

%% 

Rsite 

-site  position  vector  in  IJK 

%% 

row 

-number  of  rows  in  position  matrix 

%% 

rot  ijk  sez 

-rotation  matrix  from  IJK  to  SEZ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


format  long  g 
WGS84Data 
global  Rad 

fid  input  =  fopen  ('radardatal.dat',  'rt'); 
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i=i; 

%  scans  in  data  from  input  file  from  radar  site  or  in  this  case  the  truth  model 
while  ~feof(fid_input) 

sitlat=fscanf(  fid_input,  '%f ,  1 ) ; 

sitlon=fscanf(fid_input,  %f  ,1); 

sitalt=fscanf(fid_input,'%f,l); 

jd=fscanf(fid_input,'%f,l); 

rho  =fscanf(fid_input,'%f ,  1); 

az  =fscanf(fid_input,'%f,l); 

el  =fscanf(fid_input,'%f,l); 

range_noise=fscanf(fid_input,'%f,l); 

az_noise=fscanf(fid_input,'%f,l); 

el_noise=fscanf(fid_input,'%f,l); 

%  %  for  use  when  elevation  and  azimuth  are  given  in  degrees,  range  in  meters,  and  time  is  in 
year,  daynumber,  hour,  min, sec 

%[sitlat,sitlon,sitalt,rho,az,el,jd,range_noise,az_noise,el_noise]=unit_converter(fid_input); 

%  finds  position  vectors  from  radar  data 

[R_ijk,local_sidereal_time]  =  position_finder(jd,  sitlon,  sitlat,sitalt,  rho,az,  el); 

%  creates  matrix  of  position  vectors  as  row  vectors 
R_ijk_matrix(i,:)  =  Rijk; 

%  creates  matrix  of  j  Lilian  days 
jdtime(i,:)  =jd; 

%  creates  matrix  of  range  values 
rho_matrix(i,:)=  rho; 

%  creates  matrix  of  azimuth  values 
az_matrix(i,:)=  az; 

%  creates  matrix  of  elevation  values 
el_matrix(i,:)=  el; 

%  creates  matrix  of  1st  values 
lst(i,:)=local_sidereal_time; 

i=i+l; 

end 

%  determines  size  of  Position  matrix 
[row, column]  =  size(Rijkmatrix); 

for  ii=l:row 

%  converts  julian  days  into  time  since  initial  observation 
time  =  (jdtime(ii,:)-jdtime(151,l))*86400; 
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%  creates  matrix  of  observation  times 
timematrix(ii,:)  =time; 

end 


%%%  Linear  Least  Squares  Estimator 
%%  state  transition  matrix 

%  given  the  equations  of  motion  in  the  truth  model 

for  w  =1:300 

t  =  timematrix(w,:); 

%  defines  sub-matrices 
subO  =  eye(3,3); 
subl  =  eye(3,3)*t; 
sub2  =  eye(3,3)*(tA2)/2; 
sub3=  zeros(3,3); 
sub4  =  eye(3,3)*(tA3)/6; 


%  defines  phi_matrix 

phi  matrix  =  [subO  subl  sub2  sub4;  sub3  subO  subl  sub2;sub3  sub3  subO  subl;  sub3  sub3  sub3  subO]; 
phi_multi_array(:,:,w)=phi_matrix ; 
end 

[row2,  column,  depth]=size(phi_multi_array); 

%%  observation  relation 


%  z  =  [R_ijk]'; 

%  z  =  G(X,t)=(subO  sub3  sub3)X 
%  G=[sub0  sub3  sub3]*X; 

%  linearization  of  G 
%  H  =(partial  G)  /  (partial  X) 

H  =  [subO  sub3  sub3  sub3  ]; 

%  T=H*phi 
for  f^=l:depth 

T  =  H*phi_multi_array(:,:,f); 
T_matrix(:,:,f)  =  T; 
end 

%  differential  correction  of  intermediates 
for  s=l:row 
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K  =  [-cos(el_matrix(s))*cos(az_matrix(s)),... 

rho_matrix(s)*cos(el_matrix(s))*sin(az_matrix(s)),... 

rho_matrix(s)*sin(el_matrix(s))*cos(az_matrix(s));... 

cos(el_matrix(s))*sin(az_matrix(s)),... 

rho_matrix(s)*cos(el_matrix(s))*cos(az_matrix(s)),... 

-rho_matrix(s)*sin(el_matrix(s))*sin(az_matrix(s)); 

sin(el_matrix(s))  0  rho_matrix(s)*cos(el_matrix(s))]; 

K_array(:,:,s)=K; 


%  rotation  matrix 

rotijksez  =  [sin(sitlat)*cos(lst(s))  sin(sitlat)*sin(lst(s))  -cos(sitlat);... 
-sin(lst(s))  cos(lst(s))  0;  cos(sitlat)*cos(lst(s))  cos(sitlat)*sin(lst(s))... 
sin(sitlat)]; 

J  =  inv(rot_ijk_sez)*K_array(:,:,s); 

J_array(:,:,s)=J; 


%  instrumental  covariance 

Q_old  =  [range_noiseA2  0  0;  0  az_noiseA2  0  ;  0  0  el_noiseA2]; 
Qnew  =  J_array(:,:,s)*Q_old*J_array(:,:,s)'; 

Q_array(:,:,s)  =  inv(Qnew); 

end 


c_variance=zeros(  12,12); 


%  sums  the  values  for  the  state  covariance 
for  v  =l:row 

P=T_matrix( : , :  ,v)'*Q_array  ( : , : ,  v)  *T_matrix( : , :  ,v) ; 
c_variance=  c_variance+P; 
end 

%  inverts  the  covariance 
co  v  a  r  i  a  n  c  c~  i  n  v  ( c_  v  a  r  i  a  n  c  c ) ; 

svector=zeros(  12,1); 

%  sums  the  values  for  the  system  estimate  at  epoch  time 
for  e=  1  :row 

S=T_matrix(:,:,e)'*Q_array(:,:,e)*R_ijk_matrix(e,:)'; 
svector=  svector+S; 
end 

%  estimate  of  the  state  vector  at  epoch 
state  vector=covariance*svector 
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Appendix  D:  Estimation  Process 
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